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PREFACE 


In  recent  years,  considerable  interest  has  developed  within  the 
Air  Force  in  the  electromagnetic  pulse  generated  by  a high-altitude 
nuclear  explosion.  Possible  effects  of  this  pulse  on  military  equip- 
ment have  been  discussed,  and  there  has  been  a rekindling  of  investi- 
gations in  the  basic  physics  and  mathematics  of  the  phenomenon.  Large 
and  expensive  computer  codes  have  been  developed  to  determine  the  fields 
that  may  result.  Many  organizations  have  been  involved  in  these 
operations . 

The  Rand  Corporation  has  continuing  efforts  in  this  area,  conducted 
under  the  Project  RAND  research  project  "System  Vulnerabilities."  The 
work  has  been  communicated  to  the  Air  Force  in  several  briefings  at  the 
Air  Force  Weapons  Laboratory,  Kir t land  Air  Force  Base,  and  the  findings 
are  being  documented  in  a series  of  reports. 

This  report  presents  an  analytic  treatment  of  the  effects  of  atmos- 
pheric scattering  on  the  current  and  ionization  produced  by  the  burst. 
Studies  of  the  problem  at  other  locations  have  involved  either  ad  hoc 
physical  assumptions  or  Monte  Carlo  techniques  requiring  very  costly 
computer  programs.  The  methods  described  here  use  analytical  procedures 
as  far  as  practicable,  and  then  employ  very  rapid  computer  programs. 


Since  this  work  is  all  of  a fundamental  character,  it  is  presented  in 
unclassified  form.  Applications  to  the  electromagnetic  pulses  produced 
by  actual  nuclear  explosions  will  be  developed  in  a classified  report. 

The  analysis  and  the  computer  program  contained  in  this  report 
should  be  useful  to  persons  working  with  the  theory  or  applications  of 
electromagnetic  pulses  from  high-altitude  nuclear  explosions.  The  re- 
port presents  the  most  comprehensive  investigation  to  date  of  the  basic 

phenomena  of  current  and  conductivity.  _^CCrSSXMVJor 
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SUMMARY 


A high-altitude  nuclear  explosion  generates  gamma  rays  which  are 
absorbed  in  the  atmosphere.  The  Compton  electrons  produced  by  the 
gamma  rays  form  a current,  which  then  radiates  an  electromagnetic 
pulse,  and  also  induce  secondary  ionization  which  determines  the 
propagation  of  the  electromagnetic  field.  The  electrons  move  under 
the  influence  of  their  deflection  by  the  earth's  magnetic  field,  lose 
energy  by  the  ionization  process,  and  suffer  Coulomb  scattering  by 
the  atoms  of  the  atmosphere.  This  report  is  devoted  to  an  analytic 
investigation  of  the  current  and  ionization  produced  by  these  com- 
bined effects. 

Although  there  have  been  satisfactory  previous  treatments  of 
magnetic  deflection  and  of  energy  loss,  the  investigations  of  scatter- 
ing have  been  fraught  with  difficulties.  Methods  to  handle  scattering 
have  included  (1)  the  "obliquity"  theory,  which  contains  ad  hoc  assump- 
tions and  uncertain  physical  meaning,  and  (2)  the  Monte  Carlo  technique. 
Although  the  latter  method  is  theoretically  sound,  it  involves  lengthy 
and  expensive  computer  programs.  Accordingly,  we  have  been  impelled 
to  develop  a theory  of  scattering  which  is  carried  as  far  as  possible 
analytically. 

The  theory  is  based  on  a distribution  function  which  represents 
the  number  of  electrons  per  unit  volume  in  a unit  velocity  interval 
at  a point  in  space  and  time.  We  then  derive  a transport  equation 


which  governs  the  evolution  in  time  of  this  distribution  under  the 
indicated  physical  processes  from  the  initial  Compton  production  dis- 
tribution. This  differential  equation  is  then  solved  analytically  by 
an  expansion  to  third  order  in  the  ratio  of  ‘-spatial  spreading  to 
velocity  spreading,  thus  providing  sufficient  accuracy  to  cover  the 
significant  duration  of  the  pulse  (10  shakes),  as  is  shown  by  detailed 
numerical  investigation.  The  current  and  ionization  are  then  calcu- 
lated as  integrals  over  the  velocity  distribution  function. 

These  integrals  are  carried  out  analytically  as  far  as  we  have 
found  possible.  A computer  program  has  been  written  to  evaluate  the 


^ HffiCEDItfQ  PAGE  BLaNK-NOT  FILMED 


r 


-vi- 

integrals  which  remain.  This  program,  included  as  an  appendix,  is 
rapid  compared  to  the  Monte  Carlo  technique. 

Some  typical  results  are  presented.  For  a gamma-ray  source  taken 
to  be  a delta  function  in  time,  the  current  per  gamma  ray  at  high 
altitudes,  which  is  essentially  independent  of  scattering  and  energy 
loss,  displays  a peak  at  about  2 shakes  for  the  selected  choice  of 
magnetic  field  and  gamma-ray  energy,  and  a slow  decay  thereafter.  As 
the  altitude  decreases,  the  peak  becomes  lower,  earlier,  and  narrower, 
all  effects  which  would  be  expected  by  the  interpretation  that  scatter- 
ing reduces  the  coherence  of  the  electron  beam.  The  magnitude  of  reduc- 
tion is  considerable,  thus  at  an  altitude  of  20  km  and  times  later  than 
2 shakes,  the  current  is  reduced  by  more  than  a factor  of  four  relative 
to  the  high-altitude  current  per  gamma  ray.  The  effects  on  ionization 
are  much  less  striking,  and  the  conductivity  per  gamma  ray  associated 
with  the  ionization  is  fairly  insensitive  to  altitude. 

The  theoretical  investigation  and  computer  program  of  this  report 
should  prove  valuable  to  all  persons  interested  in  the  high-altitude 
electromagnetic  pulse. 
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I.  INTRODUCTION 


The  theory  of  the  generation  of  electromagnetic  (EM)  pulses  by 
high-altitude  nuclear  explosions  was  first  put  on  a firm  footing  in 
investigations  by  Karzas  and  Latter. ^ In  their  theory,  the  prompt 
gamma  rays  from  the  burst  interact  with  the  atmosphere  via  the  Compton 
effect  to  produce  high-velocity  electrons.  These  Compton  electrons 
are  accelerated  by  the  earth's  magnetic  field  to  produce  a radiating 
current  and  lose  energy  by  ionization  to  cause  atmospheric  conduc- 
tivity. The  propagation  of  the  fields  through  the  conducting  medium 
is  treated  by  the  "higli-f requency"  approximation,  which  neglects  cer- 
tain spatial  derivatives  relative  to  time  derivatives  and  is  regarded 
as  sufficiently  accurate  for  the  effective  duration  of  the  pulse. 

Following  this  initial  theoretical  investigation,  considerable 

effort  was  expended  in  developing  computer  codes  to  apply  the  theory 

(2) 

to  practical  cases.  Among  these  codes  were  LUMP  I,  developed  at 
the  Mission  Research  Corporation,  and  HEMP , ^ ^ HF.MP-B,^^  and  CHEMP,^^ 
developed  at  the  Air  Force  Weapons  Laboratory.  The  codes  were  gener- 
ally of  such  character  that  additional  theoretical  developments  could 
be  incorporated  without  changing  the  basic  structure  of  the  code.  As 
the  physics  of  the  calculations  have  been  improved,  the  codes  steadily 
have  become  more  complex  and  costly  to  operate. 

There  have  been  numerous  theoretical  improvements  to  the  basic 

studies.  The  dependence  of  atmnsphor \ r condor t i vi tv  on  field  strength 

has  been  treated  by  Baun/^  and  by  Higgins  et  al.^^  The  CHEMP  code^^ 

was  designed  to  include  self-consistent  electron  motions,  allowing  the 

electromagnetic  field  to  interact  with  the  electrons  which  produce  it. 

The  lag  in  the  secondary  ionization  process  was  investigated  by  Longmire 

(8) 

and  Longley.  There  have  been  several  studies  of  how  the  current 

is  affected  by  atmospheric  scattering. 

When  atmospheric  scattering  is  not  included  in  the  theory,  the 
motion  of  the  electrons  is  determined  by  their  deflection  by  the  mag- 
netic f ield  ' (electromagnet ic  field  for  self-consistent  calculations) 
and  by  the  energy  loss  to  ionization.  Thus,  there  is  a unique  and 


invertible  relation  between  the  position  and  velocity  of  an  electron 
at  a given  time  and  its  position  and  velocity  at  its  birth.  In  the 
presence  of  scattering,  an  electron  suffers  many  random  deflections 
in  the  course  of  its  motion.  Consequently,  neither  its  position  nor 
its  velocity  can  be  explicitly  deduced  from  initial  conditions.  At 
any  space-time  lo<_..  "ion  there  will  be  a distribution  of  electrons  with 
respect  to  veloci  :nd  this  enormously  complicates  the  situation. 

The  treatment  scattering  to  date  have  proceeded  along  two 
essentially  independent  lines.  First,  there  is  the  "obliquity"  theory, 
developed  in  Ref.  8 and  later  improved.  In  this  model,  the  electron 
velocity  distribution  is  replaced  by  a moving  mean  position  whose 
velocity  depends  on  the  angular  speed  of  the  motion  through  an  "obliquity" 
factor  n = < 1/cos  0 >,  where  0 is  the  angle  of  the  velocity  from  its 
initial  direction.  The  assumptions  involved  in  this  theory  have  remained 
obscure,  and  it  has  generally  been  justified  by  showing  its  agreement, 
if  any,  with  the  Monte  Carlo  theory. 

Since  the  scattering  phenomenon  is  a sequence  of  random  processes, 
it  is  amenable  to  the  Monte  Carlo  technique.  In  this  method,  the 
motion  of  each  electron  is  followed  individually.  The  scattering  is 
included  by  selecting  random  numbers  from  distributions  which  repre- 
sent the  number  of  scatterings  per  unit  path  length  and  the  angular 
deflection  suffered  at  each  scattering.  The  investigation  via  com- 
puter code  was  developed  by  Knutson^^  and  improved  by  Wittwer  et  al.^'*^ 
This  procedure  allows  a fully  satisfactory  treatment  of  scattering 
effects.  However,  it  is  extremely  costly  in  computer  time,  and  the 
change  of  any  parameter  requires  a complete  Monte  Carlo  investigation. 

The  difficulties  associated  with  these  treatments  of  scattering 
have  impelled  us  to  develop  a new  theory.  Atmospheric  scattering  is 
basically  a statistical  process.  The  obliquity  model  employs  "lumped" 
statistics,  in  that  the  characteristics  of  the  velocity  distribution 
are  lumped  into  essentially  a single  parameter,  whereas  the  Monte  Carlo 

* 

C.  L.  Longmire,  description  of  new  obliquity  given  at  a meeting 
at  Air  Force  Weapons  Laboratory,  November  1973. 
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technique  uses  individual  statistics,  treating  the  electrons  as  inde- 
pendent. We  shall  treat  the  problem  as  a general  stochastic  process, 
employing  the  methods  reviewed  and  summarized  many  years  ago  by 
Chandrasekhar . ^ ^ 

The  body  of  this  report  is  devoted  to  an  exposition  of  this  method 
of  treating  atmospheric  scattering  effects.  The  theory  is  based  on  a 
distribution  function  F(r,v,t)  which  represents  the  number  of  Compton 
electrons  per  unit  volume  per  unit  velocity  interval  at  the  position 
r at  the  time  t with  velocity  vector  v.  In  Section  II,  we  derive  the 
partial  differential  equation  which  governs  the  evolution  in  time  of 
the  distribution  function  F as  affected  by  the  various  physical  ef- 
fects. This  equation  is  referred  to  variously  as  the  Boltzmann,  Fokker- 
Planck,  or  transport  equation.  We  shall  use  the  last  term.  Effects 
included  are  production  of  electrons  by  gamma  rays,  with  the  appropriate 
Klein-Nishina  distributions  of  number  and  energy  versus  initial  direc- 
tion, electron  transport,  deflection  by  the  earth's  magnetic  field, 
energy  loss  by  ionization,  and  atmospheric  scattering.  Self-consistent 
field  effects  are  not  included,  since  they  make  the  transport  equation 
nonlinear  and  essentially  intractable.  The  current  and  the  density 
of  secondary  (ionization)  electrons  will  be  represented  in  terms  of 
the  distribution  function  F. 

Our  purpose  is  to  carry  out  the  solution  of  the  transport  equation 
by  analytic  procedures  as  far  as  possible.  These  procedures  are  de- 
veloped in  Section  III.  We  introduce  a Green's  function,  which  repre- 
sents the  distribution  of  electrons  from  an  instantaneous  uni-velocity 
point  source.  The  function  F is  then  found  as  a suitable  convolution 
of  the  Green's  function  and  the  Compton  electron  production  function. 

The  transport  equation  for  the  Green's  function  is  transformed  by  in- 
troducing new  variables  which  eliminate  all  first  derivatives  from 
the  equation.  The  resulting  equation  is  then  treated  by  successive 
approximations  in  which  the  expansion  parameter  is  essentially  the 
ratio  of  the  effect  of  positional  derivatives  to  the  effect  of  velocity 
derivatives.  This  procedure  is  carried  to  third  order.  Explicit  post 
hoc  calculations  in  Section  V show  that  the  expansion  is  sufficiently 
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accurate  for  times  up  to  10  shakes  (1  shake  = 10  ' sec),  which  covers 
the  significant  duration  of  the  pulse. 

In  Section  IV,  the  current  and  the  ionization  density  are  calcu- 
lated from  these  representations  of  the  distribution  function.  As 
many  integrals  as  possible  are  calculated  analytically.  The  remaining 
integrations,  which  only  involve  the  direction  of  the  initial  velocity, 
are  brought  into  a form  suitable  for  computer  programming. 

A computer  program  written  to  evaluate  these  integrals  proved  to 
be  very  rapid  compared  with  the  Monte  Carlo  technique.  The  program  is 
described  in  Section  V and  some  typical  results  given  and  discussed. 

The  program  itself  is  included  as  an  appendix. 

In  Section  VI,  results  are  shown  for  a typical  magnetic  field 
(0.6  gauss,  perpendicular  to  the  direction  of  the  gamma  rays)  and  a 
selected  gamma-ray  energy  (1.6  MeV) . At  high  altitudes,  the  effects 
of  scattering  and  energy  loss  are  small,  and  the  current  reduces  to  the 
Karzas-Latter  result.  It  may  be  expected  that  the  effect  of  scattering 
will  be  to  reduce  the  coherence  of  the  electron  distribution.  Accord- 
ingly, as  the  altitude  decreases,  the  peak  value  of  the  current  should 
decrease,  and  the  peak  should  occur  at  earlier  times,  since  the  electron 
beam  cannot  hold  together.  The  results  show  precisely  such  behavior. 

For  the  particular  case  shown  in  Fig.  12,  with  a delta  function  time 
dependance  for  the  gamma-ray  source,  the  current  per  gamma  ray  at  high 


altitudes  peaks  at  2.2  shakes  with  a peak  value  of  1.27  in  the  units 
of  the  figure.  At  30  km,  the  peak  occurs  at  1.6  shakes  with  a value 
of  1.06;  at  24  km  the  current  peaks  at  1.2  shakes  with  a value  of  0.69; 
and  at  18  km  it  peaks  at  0.6  shake  and  is  down  to  0.31.  Since  the  main 
electron  deposition  by  the  gamma  rays  takes  place  at  these  lower  altitudes, 
we  see  that  the  effective  radiating  source  is  very  considerably  reduced 
by  scattering,  which  will  cause  an  attendant  reduction  in  the  electro- 
magnetic field  produced  by  this  source.  Although  the  secondary  electron 
density  is  also  affected  by  scattering,  the  results  are  much  less  striking 
than  those  for  the  current.  The  steady  reduction  of  ionization  with  in- 
creasing altitude  is  primarily  because  the  energy  loss  is  proportional 
to  the  atmospheric  density.  The  conductivity  per  gamma  ray,  which  is 


I 
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proportional  to  the  number  of  secondary  electrons  divided  by  the  atmos- 
pheric density,  is  relatively  insensitive  to  altitude. 

The  calculations  of  the  radiated  field  from  the  current  and  ioniza- 
tion are  considerably  complicated  by  the  nonlinear  dependence  of  the 
effective  electron  collision  frequency  on  the  electric  field.  At  low 
altitudes,  the  effect  is  manifested  by  making  the  conductivity  field- 
dependent,  whereas  at  high  altitudes  (>  50  km),  the  effective  conduc- 
tion current  must  be  treated  by  swarm  theory.  Since  our  main  purpose 
in  this  report  is  to  present  the  effects  of  atmospheric  scattering  on 
the  collective  Compton  electron  distribution,  we  have  avoided  this 
thicket  by  not  calculating  the  fields  produced  by  the  current.  This 
makes  it  very  difficult  to  compare  our  results  directly  with  other 
investigations,  which  usually  give  only  the  radiated  fields. 

This  report  presents  our  unclassified  contribution  to  the  theory 
of  the  effects  of  atmospheric  scattering  on  the  current  and  ionization 
density  associated  with  high-altitude  nuclear  explosions.  We  have  in- 
vestigated other  phenomena  and  will  report  on  them  separately.  Applica- 
tions of  the  theory  to  the  pulses  generated  by  specific  nuclear  weapons 
will  be  presented  in  a classified  report. 
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II.  DERIVATION  OF  THE  TRANSPORT  EQUATION 

GENERAL  CONSIDERATIONS 

The  geometrical  configuration  and  coordinate  system  are  shown  in 
Fig.  1.  A nuclear  explosion  occurs  at  height  h„,  which  will  be  assumed 

D 

low  enough  that  the  curvature  of  the  earth  can  be  safely  neglected,  but 
high  enough  that  it  is  above  the  principal  altitudes  at  which  gamma 
rays  are  absorbed.  Thus,  burst  altitudes  of  100  to  perhaps  500  km  are 
considered.  This  neglect  of  the  earth's  curvature  simplifies  the  mathe- 
matics without  losing  any  of  the  essential  physics.  An  observer  is 
on  the  ground,  with  the  line  of  sight  from  burst  to  observer  making  an 
angle  ¥ from  the  vertical.  This  observer  could  also  be  in  an  airplane. 

O 


Fig.  1 — Geometry  and  coordinate  system 


A rectangular  coordinate  system  is  shown,  with  the  z-axis  along 
the  line  of  sight  directed  from  the  burst  to  the  observer.  This  will 
also  be  the  radial  direction  of  a spherical  coordinate  system.  Since 
the  effect  of  scattering  will  be  to  reduce  the  electromagnetic  pulse, 
we  shall  consider  circumstances  which  tend  to  both  simplify  the  analysis 
and  maximize  the  pulse  in  the  absence  of  scattering.  For  this  purpose. 
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we  shall  assume  that  the  earth's  magnetic  field  Is  perpendicular  to  the 
line  of  sight  and  will  define  its  direction  as  the  y-direction.  Since 
the  gamma  rays  produce  electrons  which  are  predominantly  directed  for- 
ward, this  assumption  maximizes  the  magnetic  deflection.  Self-consistent 
field  effects  are  neglected,  and  the  earth's  magnetic  field  is  assumed 
uniform  over  the  region  in  which  significant  interaction  occurs  between 
gamma  rays  and  the  atmosphere.  This  region  has  only  a 30  to  40  km 
height  range,  so  the  variation  of  magnetic  field  is  legitimately  small. 
The  methods  of  this  report  can  be  used  for  an  arbitrary  direction  of  the 
magnetic  field,  but  the  analysis  becomes  very  much  more  complicated  in 
detail.  The  coordinate  system  is  completed  by  an  x-coordinate . 

Consider  an  electron  which  is  going  to  contribute  to  the  electro- 
magnetic pulse  within  the  first  10  shakes.  This  electron  must  be  suf- 
ficiently close  to  the  line  of  sight  that  the  sum  of  the  travel  time  of 
the  gamma  ray  to  the  birth  point  of  the  electron,  the  age  of  the  elec- 
tron when  it  emits  radiation,  and  the  travel  time  of  the  radiation  to 
the  observer  does  not  exceed  by  more  than  10  shakes  the  travel  time  of 
light  along  the  line  of  sight.  As  a typical  example,  for  a burst  at 
100  km  and  an  angle  from  the  vertical  of  60  deg,  an  electron  born  at  a 
height  of  30  km  must  be  within  1.6  km  of  the  radial  axis.  Introduce 
spherical  coordinates  r.'f.X,  with  the  origin  at  the  burst  and  the  radial 
direction  along  the  line  of  sight.  Then  the  polar  angle  f will  be  less 
than  1 deg.  The  differential  equation  we  shall  derive  later  will  involve 
derivatives  with  respect  to  the  spherical  angles.  Since  the  solution 
is  nonsingular  with  respect  to  angles  and  the  total  angular  variation 
is  small,  we  shall  neglect  these  spherical  angle  derivatives  of  the 
electron  number  density  compared  to  the  radial  derivative.  We  have 
worked  out  the  extremely  tedious  calculations  required  by  not  making 
this  approximation  and  find  that  the  current  and  conductivity  are  only 
modified  by  a few  percent  over  the  time  interval.  Furthermore,  the 
distance  traveled  by  an  electron  in  the  10  shake  retarded  time  interval 
is  small  compared  to  the  atmospheric  scale  height,  so  we  shall  neglect 
the  variation  of  the  atmospheric  density  over  this  distance.  We  thus 

neglect  effects  on  the  order  of  ct  / H,  where  c is  the  velocity  of  light, 

m 

T the  pulse  duration,  and  H the  scale  height, 
m 
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The  gamma  rays  produce  electrons  In  all  directions.  Hence,  the 

-f 

velocity  vector  v shown  in  Fig.  1 can  have  all  possible  directions 
and  magnitudes  up  to  the  maximum  permitted  by  the  Compton  effect.  The 
rectangular  components  of  the  velocity  (v^,v  ,vz)  are  given  in  terms 
of  its  magnitude  v = 6c  and  direction  angles  0,(p,  by 


vx  = 6c  sin  0 cos  <p 


v = 6c  sin  0 sin  cp 

y 


v = 0c  cos  0 
z 


(2.1a) 

(2.1b) 

(2.1c) 


The  polar  angle  0 is  indicated  in  Fig.  1. 

We  shall  now  introduce  stochastic  procedures,  following  the  method 
described  on  p.  42  of  Ref.  11.  An  electron  moves  under  the  action  of 
the  coherent  acceleration  processes  (i.e.,  energy  loss  and  magnetic 
deflection)  and  the  incoherent  scattering  process.  We  assume  a time 
interval  At  exists  such  that  the  electron  undergoes  many  infinitesimal 
scatterings,  but  that  its  position  and  velocity  do  not  change  appreciably. 
Then  the  increments  in  velocity  and  position  in  this  time  interval  are 
given  by: 


Av  = aAt  + 6v 

Ar  = vAt 


(2.2a) 

(2.2b) 


where  a is  the  acceleration  produced  by  the  coherent  processes,  and 
6v  is  a fluctuating  quantity  with  a definite  law  of  distribution.  Thus, 
the  transition  probability  that  an  electron  with  velocity  v at  position 
r will  undergo  a velocity  change  Av  and  position  change  Ar  is  given  by 
the  function: 


¥(r,v,Ar,Av)  - f(r,v,Av)  6(Ar  - vAt) 


(2.3) 


Here,  6 denotes  the  Dirac  delta  function  and  is  not  to  be  confused  with 
the  fluctuation  velocity  6v. 
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Let  us  now  define  the  distribution  function  F such  that 

F(r,v,t)  drdv 


is  the  expected  number  of  electrons  at  time  t in  an  element  of  spatial 
volume  dr  in  the  vicinity  of  r and  in  an  element  of  velocity  volume  dv 
in  the  vicinity  of  v.  We  ask  how  many  electrons  are  in  this  phase  space 
volume  at  time  t + At . Clearly,  this  is  the  sum  of  those  produced  by 
the  gamma  rays  in  the  volume  during  the  time  At  and  the  net  change  pro- 
duced by  inf low/ outflow  from  the  volume  as  a result  of  the  acceleration 
forces.  Thus  (see  Ref.  11,  Eq.  (3.24)),  we  have  the  integral  equation: 


F(r,v,t+At)  drdv  = Ne(r,v,t)  Atdrdv 


//' 


F(r-Ar , v-Av, t)  ^(r-Ar ,v-Av,Ar ,Av)  dArdAvdrdv 


(2.4) 


Here,  is  the  rate  of  production  of  electrons  by  gamma  rays,  and  the 
integral  term  depicts  transitions  in  and  out  of  the  phase  space  volume. 
The  representation  (2.3)  permits  integration  over  dAr.  With  a slight 
change  of  variables,  Eq.  (2.4)  takes  the  form: 


F(r+vAt , v+aAt , t+At)  drdv  = N (r,v,t)  Atdrdv 


/' 


+ / F(r,v-6v,t)  'Kr ,v-6v,Av)  dAvdrdv 


(2.5) 


Under  the  previous  assumptions,  the  spatial  dependence  of  F involves 
the  radial  coordinate  only.  The  functions  on  the  left  and  right  sides 
of  Eq.  (2.5)  can  be  expanded  in  Taylor  series,  and  the  limit  At  -*■  0 
taken.  The  transport  equation  results,  which  governs  the  evolution 
in  time  of  the  distribution  function  F: 
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ff-  + v cos  0 -2-  1^-  r2F  + V+  • aF  = Nfi(r,v,t) 

r it  /*2u 

+ Nv  | I sin  Q'  d0  'dcp'Q(B,0,cp,0",cp')  [F(r,B,0',<p',t) 

J o J o 

- F (r , B,  0,cp,  t)  ] (2.6) 

Here,V^  is  the  gradient  operator  in  velocity  space,  N is  the  number  of 
atoms  per  unit  volume,  and  Q is  the  differential  cross  section  for 
elastic  scattering  from  a velocity  described  by  the  direction  angles 
Q' ,<p'  to  the  angles  0,<p.  The  transition  probability  T is  the  product 
NvQ.  The  two  terms  in  brackets  in  the  integral  indicate  that  scatter- 
ing takes  place  both  into  and  out  of  the  element  of  volume  in  velocity 
space.  The  space  divergence  operator  on  the  left  represents  the  net 

flow  through  the  spatial  element  of  volume,  taking  into  account  the 

2 

change  in  the  element  of  surface  area  r du),  where  du)  is  an  element  of 
solid  angle. 

Equation  (2.6)  is  the  fundamental  transport  equation.  We  shall 
now  consider  the  explicit  representation  of  the  coherent  acceleration, 
scattering,  and  production  terms. 

t 

COHERENT  ACCELERATION 

The  coherent  acceleration  is  given  by  the  magnetic  deflection  and 
the  energy  loss  expressions.  The  equation  of  motion  of  a relativistic 
electron  under  these  forces  is: 


d mv 

dt  (l-B2)*5 


ev  x B - W(8) 


v 

V 


(2.7) 


Here,  B is  the  velocity  of  the  electron  divided  by  c,  the  velocity  of 
light,  B is  the  earth's  magnetic  field,  and  W(B)  is  related  to  the 
energy  loss  in  a manner  to  be  demonstrated.  Self-consistent  field 
effects  have  been  omitted  in  Eq.  (2.7).  If  we  take  the  scalar  product 
of  v and  Eq.  (2.7)  and  simplify,  there  results: 


d_ 

dt 


(2.8) 


" 
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me 


(1-B2)*4 


vW(B) 


The  left  side  is  the  rate  of  change  of  the  electron  energy.  Conse- 
quently, the  function  W(0)  is  given  by 


W(B) 


dE 

dS 


(2.9) 


Thus,  W(B)  is  the  energy  loss  per  unit  path  length  dS.  This  is  the 
quantity  calculated  and  tabulated  in  the  literature. 

We  have  assumed  here  that  all  the  inelastic  scattering  processes 
can  be  represented  as  an  energy  loss  term  which  is  parallel  to  the 
electron  velocity.  Actually,  these  inelastic  processes  also  involve 
change  in  direction,  but  such  effects  are  quite  small.  The  dominant 
energy  loss  process  at  Compton  energies  is  ionization,  with  small 
momentum  change  per  collision.  The  change  in  the  direction  of  the 
velocity  of  the  primary  for  such  ionizing  collisions  has  been  shown 
experimentally  to  be  small.  We  neglect  the  very  rare  large-angle  in- 
elastic scattering  events.  With  these  assumptions,  the  use  of  the 
term  "coherent"  for  energy  loss  seems  satisfactory.  This  treatment  is 
consistent  with  all  other  investigators. 

When  Eq.  (2.8)  is  used  in  (2.7),  we  obtain  for  the  coherent  accel- 
— ► 

eration  a the  form: 


S - - S (l-B2)*5  v * s - «<b)»-b2)3/2; 

m ^ « 

me  g 


(2.10) 


This  expression  is  to  be  substituted  into  Eq.  (2.6).  The  first 
term  is  perpendicular  to  the  velocity,  and  the  divergence  will  only 
involve  derivatives  with  respect  to  the  angles  0,cp.  The  second  term 
is  parallel  to  the  velocity,  and  only  derivatives  with  respect  to  the 
velocity  magnitude  will  appear.  We  have  assumed  the  magnetic  field 
to  be  in  the  y-direction.  When  the  derivatives  are  simplified,  we  have: 
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aF  = 


- — (1-B2) 

m 


\ T 3F  , cos  0 sin  w 9F 

- cos  cp  — + — 


90 


sin  9 3cp. 


_ l-i-f 
a2  9B  L 


62(1_62)3/2  W(B)  p 
me 


(2.11) 


where  B is  the  magnitude  of  the  earth's  field  and  we  have  replaced  the 
derivative  with  respect  to  v by  that  with  respect  to  B. 

The  function  W(8),  the  energy  loss  per  unit  path  length,  has  been 

( i ) ) 

calculated  by  Bet he;  “ it  is  given  by: 


M(B) 

me 


2„NZr2c  M<“>2 

e_  ’ 

3 - (~  “ ~~2>  l°g  2 + ~2  + 8 1 

T Y Y 

Y = (I-bV^2 


(2.12) 


(2.13) 


Here,  N is  the  number  of  atoms  per  unit  volume,  Z is  the  number  of 
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electrons  per  atom  (Z  = 7.24  for  air),  and  r^  = e /me  . The  quantity 

I is  an  effective  ionization  potential  which  is  80.5  eV  for  air.  The 

2 2 
total  energy  of  the  electron  is  ymc  and  its  kinetic  energy  is  (Y~l)  me  • 

The  expression  (2.12)  has  dimensions  of  inverse  time.  We  shall 

measure  time  in  shakes,  so  the  multiplying  factor  becomes 

2TTNZr^c  = 0.05522  p [shakes-1]  (2.14) 

The  quantity  p = P/PQ  is  the  ratio  of  the  air  density  at  the 
altitude  at  which  the  electrons  are  located  to  the  density  at  sea 
level.  The  factor  in  brackets  in  Eq.  (2.12)  will  be  designated  V(8) 
and  is  plotted  against  electron  kinetic  energy  in  Fig.  2,  where  we 
also  plot  V(B)/8.  It  may  be  seen  from  Fig.  2 that  V(B)/B  only  varies 
slightly  over  most  of  the  range  in  energy.  If  it  is  replaced  by  the 
constant  value  21,  the  error  will  not  exceed  3.5  percent  for  energies 
between  250  KeV  and  2 MeV,  and  will  only  be  20  percent  at  100  MeV. 
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Electron  kinetic  energy,  MeV 


Fig.  2 — Energy  loss  factor  versus  energy 
Thus,  we  shall  use  the  approximate  form: 

W(6)^  = 1 -16-P  [ shakes'1 1 (2.15) 

me  p 

The  coefficient  of  the  magnetic  field  term  is 

— = 0.17588  B (gauss)  [shakes  1]  (2.16) 

m 

The  form  (2.15)  for  the  energy  loss  per  unit  path  length  is  equiva- 
lent to  a constant  loss  per  unit  time,  so  the  energy  decreases  linearly 
with  time  in  the  reference  frame  of  the  electron,  down  to  relatively  low 


-14- 


energies.  At  these  low  energies,  the  electron  loses  energy  much 
Caster  and  essentially  stops.  We  shall  use  the  form  (2.15)  for  all 
energies  and  will  not  expect  to  affect  the  current  or  ionization 
significantly. 


ATMOSPHERIC  SCATTERING 

The  atmospheric  scattering  term  on  the  right  hand  side  of  Eq.  (2.6) 
involves  the  differential  cross  section  for  elastic  scattering.  This 
will  generally  be  a function  of  the  change  in  the  momentum  of  the  elec- 
tron. As  discussed  above,  we  shall  neglect  all  inelastic  scattering 
processes,  which  will  be  regarded  as  represented  by  the  energy  loss. 

The  change  in  momentum  involves  the  invariant  magnitude  of  the  velocity, 
and  the  cosine  of  the  angle  between  the  initial  and  final  direction  of 
the  velocity.  We  thus  have: 


Q(B,0,cp,0^,cp*)  = Q[B,  cos  0 cos  Q'  + sin  0 sin  0*  cos  (tp-cp*)  ] (2.17) 


We  shall  call  the  second  argument  cos  x-  The  function  Q is  strongly 
peaked  at  small  values  of  x*  We  therefore  expand  the  function 
F (r , 0, 0 ' ,cp' , t ) in  a Taylor  series  in  Q'-O,  cp^-c p.  To  second  order,  we 
have : 


F(r,0,0',cp',t)  = F(r,0,0,<p,t)  + ^ (0'-0)  + g—  (cp'-cp) 


(2.18) 


2 

1 3 F 


, a o r „s  2 , 9 F ax  , - x , 1 9 F y > x 2 

+ 2 ^2  (9  “0)  + 909^  (0  "e)(cP  ^ + 2 ^2  (<P 


Since  the  integration  is  over  the  direction  cosines  associated 
with  0',cp',  the  angle  differences  should  be  represented  in  terms  of 
differences  of  direction  cosines.  Again  to  second  order: 


0-0 


„ _ (cos  0 - cos  Q')  cos  0 (cos  9"  - cos  9)' 


iin  9 


2 sin  0 


(2.19a) 


sin  Q'  sin  (cp'-<p) 
sin  9 


1 + 


sin  9 - sin  d'  cos  ( cp'-cp ) 
sin  0 


— 


cp  -t p 


(2.19b) 
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These  expressions  are  to  be  inserted  into  Eq . (2.18),  then  (2.18) 

itself  is  to  be  put  into  the  integral  of  Eq . (2.6).  With  cp'-cp  as  a 

2 

new  azimuthal  angle,  the  terms  in  9F/09  and  9 F/909cp  are  odd  functions 
and  integrate  to  zero.  The  integration  over  the  direction  of  v can  be 
simplified  by  rotating  the  coordinate  system  so  the  direction  0,tp  repre- 
sents the  polar  direction.  In  this  coordinate  system,  the  angle  x» 
defined  following  Eq . (2.17),  is  the  polar  angle.  We  will  call  the 
corresponding  azimuthal  angle  a.  The  relation  between  the  quantities 
required  in  Eq.  (2.19)  for  the  0',tp'  and  x>a  systems  is: 

cos  0 - cos  Q'  = cos  0 (1  - cos  x)  + sin  9 sin  X cos  a (2.20a) 

sin  0 - sin  0'  cos  (cp'-cp)  = sin  0 (1  - cos  x)  - cos  0 sin  X cos  a (2.20b) 

sin  0"  sin  (cp"-cp)  = sin  X sin  a (2.20c) 

Since  Q does  not  involve  a,  the  integration  over  a is  immediate. 

To  this  point,  the  scattering  integral  of  Eq.  (2.6)  becomes: 


! 
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We  carried  this  procedure  to  sixth  order;  the  result  is: 


S2  + 32  S4  + k h)  V'F 


+ {h s*  * nfe  SJ  ^ 2p  * ms  h vV,2p 


sin*"11^1  x Q(B,  cos  x)  dX 


(2.23a) 


(2.23b) 


The  moments  S cun  be  expected  to  fall  rapidly  with  m.  We  now 
require  .in  expression  for  t lie  elastic  scattering  cross  section  Q.  We 
are  concerned  with  electron  energies  such  that  the  momentum  transfer 
is  smn  1 1 m.i  the  :i.  rn  approximation  to  scattering  is  valid.  From  the 
review  irt  i.  li  \ Mel.  et  al.,  we  have  selected  the  formula  of 

M l i,  re,  ."ii.  •.•pi,  , , ni  - in  improved  version  of  the  Fermi-Thomas 

theorv.  This  i rmula 


AzV 


He,  CO*  |) 


' Et- 

-8  *7*  a:  + 


1 i 


(2.24a) 


,1/3 


sin  \/2 

(2.24b) 

b / 121 
i 

(2.24c) 

0.55  a3  = 0.35 

(2.24d) 

1.2  b3  = 0.3 

(2.24e) 

This  expression  may  be  substituted  into  Eq.  (2.23b)  for  the 
scattering  moments.  All  integrals  may  be  performed  analytically.  The 
result  is  an  expansion  in  A^/28y*  beginning  with  a logarithmic  term 
(in  only),  then  constants,  then  positive  even  powers.  This  expan- 

sion parameter  will  be  small  for  all  significant  values  of  8,  so  only 
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the  logarithm  and  constants  will  be  kept.  To  this  order: 


s _ 4z2r2  Ozf> 

2 e g4 


!°g  ^f/31  " 1-0793 


S = 4Z2r2  -(1  ^ ) 2/3 
4 e 64 


, r.2  .2  (1-3  ) 

S.  = 4Z  r 2/5 

6 e 04 


(2.25a) 


(2.25b) 


(2.25c) 


3 


We  shall  modify  these  formulas  by  replacing  Z by  Z(Z+1).  This 

is  the  "standard"  procedure  for  including  scattering  by  the  orbital 

electrons,  whereas  Eq.  (2.24)  corresponds  to  scattering  by  the  screened 

nucleus.  When  Eqs.  (2.25)  are  substituted  into  (2.23a),  the  coeffi- 
2 

cient  of  V F becomes 


2 

2TTZ(Z+l)r2Nc  log  (114.78  ByZ“1/3) 

6 3 


(2.26) 


If  the  higher  order  terms  had  not  been  included,  the  coefficient  in 
the  argument  of  the  logarithm  would  have  been  82.24.  The  logarithm  is 
between  3 and  5 over  the  energy  range.  The  coefficients  of  the  higher 
derivatives  in  Eq.  (2.23a),  corresponding  to  the  logarithm  in  Eq.  (2.26), 
are  respectively  0.05972  and  0.00069,  so  it  is  fully  legitimate  to 
neglect  the  fourth  and  sixth  derivative  terms.  Substituting  numbers 
into  Eq . (2.26)  yields: 


0.455  p log  (59.33  3y) 

3 


[shakes  3] 


(2.27) 


See  Ref.  12,  pp.  261  and  284  for  many  other  references. 
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The  logarithm  and  log/B  are  plotted  against  energy  in  Fig.  3. 
Again,  log/B  varies  only  slightly  over  the  entire  range.  We  shall 
replace  it  by  the  constant  value  5.52,  which  is  in  error  by  only  3 
percent  between  250  KeV  and  2 MeV,  and  only  20  percent  at  100  KeV. 

The  ratio  between  the  scattering  coefficient  and  the  energy  loss  co- 
efficient is  more  nearly  constant  than  either  of  them.  This  was 
initially  pointed  out  in  Ref.  8.  (The  argument  of  the  log  was  65.36 
By  in  Ref.  8,  representing  a slightly  different  choice  of  lower  limit 
in  the  integration.)  With  this  constant  value,  the  scattering  inte- 
gral reduces  to: 


2.51  p V2F  [shakes'1]  (2.28) 

B 

This  will  be  used  in  Eq.  (2.6)  at  all  energies. 


Fig.  3 — Scattering  coefficient  versus  energy 


-19- 


COMPTON  ELECTRON  PRODUCTION 

The  term  in  Eq.  (2.6)  represents  the  production  of  Compton 
electrons  by  gamma  rays  and  acts  as  a driving  function  for  the  differ- 
ential equation.  The  gamma  rays  will  interact  with  the  atmosphere  via 
the  Compton  effect,  and  the  number  of  electrons  produced  per  unit  time 
will  be  proportional  to  the  number  of  gamma  rays  which  undergo  inter- 
action and  to  the  differential  cross  section  for  production  of  an 
electron  with  a specified  momentum  vector.  At  each  scattering  of  a 
primary  gamma  ray,  a secondary  gamma  ray  will  also  be  produced,  which 
will  in  turn  be  scattered  at  a later  time,  leading  to  a general  cascade 
process.  The  retarded  time  until  a second  scattering  takes  place  will 
usually  exceed  the  10-shake  retarded  time  scale  of  this  report,  so  we 
will  limit  our  considerations  here  to  the  production  of  electrons  by 
the  first  scatterings  of  gamma  rays  from  the  burst.  We  have  investi- 
gated the  cascade  process  on  a longer  time  scale  and  will  report  the 
results  in  a separate  publication. 

We  shall  consider  a source  of  gamma  rays  of  a single  energy, 

2 

E = kmc  . If  the  burst  is  represented  by  an  energy  distribution, 

our  results  can  be  integrated  over  k to  obtain  the  total  current  and 

ionization.  The  nature  of  the  Compton  scattering  process  is  such  that 

there  is  a unique  relation  between  the  angle  at  which  an  electron  is 

emitted  and  its  velocity.  The  function  N can  then  be  written  as  a 

e 

product  of  four  factors: 


N (r  v,t)  = f(t  ~ r/c)  K(cos  6)  6 (B  - h(cos  6)) 

6 ’ ’ 2tt$2 


(2.29) 


These  functions  have  the  following  interpretations.  First,  g(r) 
is  the  total  number  of  gamma  rays  which  are  absorbed  per  unit  volume 

at  the  distance  r from  the  burst.  Let  0 be  the  total  Compton  scatter- 

(15)  C 

ing  cross  section,  given  by: 


2 2 + 8k  + 9k2  + k3  _ (2  + 2k  - k2)  lot; 

c 6 k2 (1  + 2k)2  2k3 


(1  + 2k) 


(2.30) 
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» 


Let  the  yield  of  the  burst  in  MeV  of  gamma  rays  be  Y.  Then  the  func- 
tion g(r)  is  given  by 


where  s(k)  is  the  function  in  brackets  in  Eq.  (2.30),  Nq  is  the  atomic 
density  at  sea  level,  and  p is  the  ratio  of  the  density  at  r to  the 
sea  level  density.  Equation  (2.31a)  involves  a mixture  of  units.  With 
some  further  changes  of  units,  we  will  write 


(2.32) 


where  now  Y^  is  the  gamma  yield  in  kilotons,  the  distance  r from  the 
burst  is  measured  in  kilometers,  the  gamma  energy  is  in  MeV,  and  the 
total  deposition  g(r)  is  per  cubic  centimeter.  The  integral  in  the 
exponent  is  proportional  to  the  total  air  mass  traversed  by  the  gamma 
rays  to  reach  the  distance  r.  We  have  plotted  g(r)  versus  altitude 
in  Fig.  4 for  a burst  at  100  km  altitude,  a yield  of  1 kT,  a gamma 
energy  of  1.6  MeV  (k  = 3.1311),  and  an  angle  from  the  vertical  of  60 
deg.  At  this  energy,  s(k)  = 0.3321.  The  production  for  other  yields, 
energies,  and  geometries  will  have  the  same  general  character,  a peak 
at  some  altitude  in  the  20  to  40  km  range  with  a sharp  drop-off  at 
lower  altitudes,  and  an  exponential  decline  at  higher  altitudes  es- 
sentially proportional  to  density.  We  have  used  the  Cospar  Interna- 
tional Reference  Atmosphere  1961 . ^ ^ The  peak  is  actually  quite  broad, 
becoming  half  the  value  at  the  height  of  maximum  (33  km)  at  25.5  and 
at  44.5  km.  The  density  varies  by  a factor  of  37  over  this  altitude 
range . 
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Altitude,  k™ 


Fig.  4 — Total  electron  production  versus  altitude 

The  function  f ( t — r/c)  represents  the  time  history  of  the  gamma 
source.  Rather  than  use  specific  functions,  we  shall  calculate  with 
a delta  function  for  f.  The  actual  time  dependence  of  the  source  can 
then  be  brought  in  by  convolution.  The  computer  program  in  the  appen- 
dix can  be  augmented  by  a simple  integration  process  to  obtain  the 
current  and  ionization  produced  by  an  arbitrary  gamma-ray  source  time 
dependence.  No  further  differential  equations  need  be  solved.  We 
anticipate  that  any  actual  time  dependence  of  the  gamma  source  will 
cause  the  current  to  display  a slower  initial  rise  rate  than  will  be 
obtained  with  the  delta  function. 

The  function  K(cos  0)  corresponds  to  the  angular  distribution  of 
the  Compton  production  process.  The  literature  generally  gives  this 
Klein-Nishina  distribution  in  terms  of  the  differential  cross  section 
for  producing  a gamma  ray  of  momentum  k'  from  a gamma  ray  of  momentum 
k.  The  cross  section  for  electron  production  is  obtained  from  this 


expression  by  replacing  1c  by  k - p,  where  p is  the  electron  momentum. 
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and  then  multiplying  by  the  ratio  of  the  unit  solid  angle  for  the 
”►  > -> 

photon  k to  that  of  the  electron  p.  We  normalize  so  the  integral  of 
K(cos  0)  sin  0 is  unity.  Omitting  the  detailed  algebraic  manipulations, 
there  results: 


K(x) 


2 (1+k) 2x 

s(k)[(l+k)2  - k2x2]2 


(1+k)2  - k(2+k)x2  + (1+k)2  - k2x2 

(1+k)2  - k2x2  (1+k)2  - k(2+k)x2 

2 2 2 

_ 4 (1+k)  x (1-x  ) 

[(1+k)2  - k(2+k)x2]2 


(2.33) 


where  x = cos  0.  We  have  plotted  K(cos  0)  and  K(cos  0)  sin  0 in  Figs. 

5 and  6 for  E = 1.6  MeV  (k  = 3.1311).  These  are  respectively  the 
number  per  unit  solid  angle  and  the  number  per  unit  angle.  Note  the 
factor  of  10  change  in  scale  between  Figs.  5 and  6.  The  sharp  peak  of 
K(cos  0)  at  small  angles  is  evident  in  Fig.  5.  However,  the  curve  of 
Fig.  6,  the  number  versus  angle  distribution,  which  will  appear  in  the 
calculation  of  current  and  ionization,  begins  at  zero,  peaks  at  10  deg, 
and  falls  off  essentially  linearly.  This  much  broader  distribution 
indicates  that  the  often-made  assumption  that  the  Klein-Nishina  dis- 
tribution can  be  replaced  by  forward-directed  electrons  is  fraught  with 
considerable  danger.  Only  half  the  number  of  electrons  are  produced  at 
angles  below  27  deg. 

The  final  delta  function  factor  in  Eq.  (2.29)  relates  the  electron 
velocity  to  angle.  The  normalization  is  such  that  the  integral  over 
B will  be  unity.  The  electron  velocity  and  kinetic  energy  are  found 
from  the  Compton  formulas,  and  are: 


h(cos  0) 


2k(l+k)  cos  0 
(1+k)2  + k2  cos2  0 


E(cos  0) 


2 oi.2  2 

me  2k  cos  0 

2 2 2 
(1+k)  - k cos  0 


(2.34) 


(2.35) 
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These  functions  are  plotted  against  angle  in  Fig.  7 for  E = 1.6  MeV. 
Although  the  velocity  drops  slowly  with  angle,  the  energy  drops  rather 
steeply.  The  average  energy  per  electron  can  be  calculated  using 
Figs.  6 and  7;  it  is  almost  exactly  0.8  MeV.  Thus,  half  the  energy 
of  the  original  gamma  rays  has  gone  into  the  production  of  these  first- 
generation  electrons  and  half  remains  in  the  photon  cascade  for  later 
deposition.  Half  the  total  energy  of  the  distribution  is  associated 
with  electrons  produced  at  angles  less  than  18  deg,  so  the  effective 
energy  distribution  E(cos  0)  K(cos  0),  not  shown,  is  narrower  than  the 
electron  number  distribution. 


Fig.  7 — Electron  velocity  and  kinetic  energy;  Ey  = 1.6  MeV 

CONNECTION  WITH  CURRENT  AND  SECONDARY  IONIZATION 

At  this  point,  we  can  write  the  complete  equation  for  the  evolu- 
tion of  the  distribution  function  F.  Using  the  results  of  the  previous 
subsections,  we  have: 


-25- 


Tt  + Be  co.  6 ij  |j  r2F  + f d-S2)*5  (cos  <p  f - 


cos  6 sin 
sin  0 


J2  9F\ 
3cp/ 


- i-  6(l-62) 3/2  F - 2.51  p(r)  ii=P-  M 


36 


30 


cos  0 3F 

• A f\rt  ' 


i rr 

2 


sin  0 30  2 . „ 

sin  0 3cp 


g(r)  f (t  - r/c)  K(cos  6)  6 ( B — h(cos  9)) 

2 

2tt(3 


(2.36) 


This  is  the  complete  transport  equation.  The  time  is  in  shakes  and, 
if  Eq.  (2.32)  is  to  be  used  for  the  total  production,  the  units  of  F 
are  electrons  per  cubic  centimeter.  The  reduced  density  p(r)  is  a 
function  of  altitude  or,  equivalently,  of  range  to  the  burst. 

The  physical  quantities  to  be  determined  from  F are  found  from 
Maxwell's  equations.  In  MKS  units,  these  are: 


V x E = - 


3B 

3t 


(2.37a) 


C 


(2.37b) 


where  J denotes  the  total  current,  including  the  Compton  current 
and  the  conduction  current  Jg.  At  altitudes  below  50  km,  the  secondary 
electrons  experience  so  many  collisions  with  the  neutral  molecules  of 
the  atmosphere  in  a time  during  which  the  electromagnetic  field  changes 

-f 

appreciably  that  the  conduction  current  is  a true  conduction  term  OE. 
The  conductivity  is  then  given  by: 


e Ns(r,t) 
mv(r) 


O 


(2.38) 
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where  v(r)  is  the  collision  frequency  which  is  proportional  to  the 
density  and  N^(r,t)  is  the  total  number  of  secondary  electrons  gen- 
erated in  the  atmosphere  by  the  ionization  process. 

The  Maxwell  equations  are  reduced  and  simplified  by  dropping  the 
angular  derivatives  and  employing  the  high-frequency  approximation  of 
Ref.  1.  This  enables  us  to  write  a single  equation 

f-  rE  + 222  rE  = 22  rJ  (2.39) 

or  2 2 c 


where  the  arguments  of  E and  are  r,  t - r/c.  The  derivative  is 
with  respect  to  r,  with  t - r/c  held  constant.  Only  the  transverse 
components  of  E or  appear  ia  Eq . (2.39).  Thus,  the  quantities  we 
shall  calculate  from  F are  yc.fc/2  and  Ng.  We  note  that,  although  v(r) 
has  been  written  as  a function  of  position  only,  it  actually  is  strongly 
dependent  on  the  field,  so  Eq . (2.39)  is  a nonlinear  equation  for  the 
field  which  must  be  integrated  step-vise.  Since  we  shall  not  calculate 
fields  in  this  report,  we  shall  not  treat  Eq.  (2.39)  further. 

The  Compton  current  is  obtained  from  the  velocity  distribution 
function  by  multiplying  by  ev  and  integrating  over  all  values  of  the 
velocity.  This  yields: 


2 c 


0 cos  9 dGdcp  F(r ,S,6,cp,t) 


(2.40) 


—8 

The  numerical  factor  e/2e  = 0.90467  * 10  . The  other  component  will 

integrate  to  zero. 

For  the  secondary  electrons,  we  have  the  rate  relation 


E^N 


S 


(2.41) 


which  asserts  that  the  number  of  secondaries  produced  per  unit  time 
multiplied  by  the  mean  energy  per  secondary  equals  the  number  of  pri- 
maries multiplied  by  the  energy  loss  per  unit  time  of  the  primaries. 
Using  a mean  energy  per  secondary  of  34  eV  and  with  the  help  of  Eq. 
(2.15),  we  have 
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<5NS  = 17435  pF 


(2.42) 


The  total  rate  of  production  of  secondaries  is  thus: 


N 


S 


17435 


sin  0 d0d<p  F(r,0,0,cp,  t) 


(2.43) 


The  integral  on  the  right  is  also  the  total  number  of  primaries  N^. 

To  obtain  the  actual  number  of  secondaries  from  their  rate  of 
production,  we  shall  use  the  model  of  Ref.  8,  which  takes  into  account 
the  manner  in  which  the  electrons  produced  by  Compton  ionization  in 
turn  produce  additional  ionization.  With  this  model,  we  have: 


Ng  = f dt'N  (t')  L(t-t') 

J —CO 

L(t-t  ')  = ||  J^l  + 0.16  log  ( 1 + 331 00  p(t-t ')) 


= 1 if  Eq.  (2.44b)  exceeds  1 


(2.44a) 

(2.44b) 

(2.44c) 


We  have  modified  the  actual  formula  of  Ref.  8 by  introducing  p in- 
stead of  the  atomic  density  N and  have  expressed  time  in  shakes.  Elec- 
tron attachment  has  been  neglected,  which  should  be  satisfactory  at 
these  altitudes  on  this  time  scale.  We  thus  see  that  the  quantities 
to  be  calculated  from  F are  the  integrals  of  Eqs.  (2.40)  and  (2.43). 

With  Eqs.  (2.36),  (2.40),  and  (2.43),  we  have  completed  the  deri- 
vation and  shall  now  proceed  to  the  solution.  Lengthy  and  complicated 
mathematical  treatment  follows.  We  considered  delegating  this  material 
to  appendices,  but  feel  the  report  would  be  Impossible  to  follow  if  we 
leaped  from  the  formulation  of  the  transport  equation  to  (1)  its  solu- 
tion, which  in  abbreviated  form  in  Eq.  (3.27)  occupies  ten  lines  of 
solid  mathematical  symbols,  or  (2)  the  representations  for  current  and 
ionization  rate,  which  are  of  comparable  length.  We  have  attempted  to 
make  the  presentation  clear  without  becoming  either  too  concise  or 
hopelessly  bogged  down  in  detail. 
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III.  SOLUTION  OF  THE  TRANSPORT  EQUATION 

THE  GREEN'S  FUNCTION 

Equation  (2.36)  represents  the  evolution  in  time  of  the  electron 
distribution  function  in  response  to  the  Compton  production  function. 

We  shall  introduce  a Green's  function,  which  is  the  response  to  an 
instantaneous  point  source  with  a specified  velocity  vector.  The  true 
distribution  function  can  then  be  determined  by  a convolution  between 
the  Green's  function  and  the  Compton  production  function.  The  equa- 
tion for  the  Green's  function  is  simpler  to  solve  than  is  Eq.  (2.36). 

First,  we  shall  consider  the  distance  scale  of  the  electron  motion. 
In  the  absence  of  scattering,  the  fastest  electrons  of  the  initial  dis- 
tribution will  spiral  around  the  field  lines  with  a helical  radius  of 
about  100  m.  We  do  not  expect  scattering  to  change  this  scale  signifi- 
cantly. This  scale  is  small  compared  to  the  6 km  atmospheric  scale 
height,  so  we  shall  neglect  the  variation  of  the  density  p(r)  and  the 
deposition  g(r)  over  the  electron  parh.  We  extract  a velocity-dependent 
factor,  and  write: 


_ g(r)  W(r ,B,9,cp,t) 

B(i-B2)3/2 


(3.1) 


9W 

at 


. n aw  . eB  ,,  / aw  cos  0 sin  cp  3W  \ 

+ Be  cos  <p  ^ + - (1-6  ) (ops  cp  ^ s-in  e V ^ ) 


1.16  p M . 2.51 5 (U|!)  ,2„ 


,2,3/2 


^ f(t  - r/c)  K(cos  0)  6 (b  - h(cos  0)) 


(3.2) 


j 

2 J 


c(r'’  / t~~2~St2  / / sln2  e' cos  * d0'dcp' 

Jo  (1-8  ) Jo  Jo 


W(r',8',d',t p',t') 


(3.3) 
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N„(r',t')  = 17435  pg(r')  f — — -d-| f f sin  0'  dO'dcp' 

Jo  (1-r V'Vo  Jo 


W(r",0',0',cp',t') 


(3.4) 


where  the  primes  have  been  introduced  to  facilitate  later  work. 

The  point  of  the  transformation  (3.1)  is,  of  course,  to  eliminate 
the  velocity  factor  under  the  6 derivative.  The  initial  condition  on 
the  distribution  W is  that  W vanish  for  t < r/c. 

The  Green's  function  G,  a function  of  the  10  variables,  r,0,0,cp,t, 
is  taken  to  satisfy  the  differential  equation 


3G 

3t 


+ 0c  cos  0 


9G 

3r 


+ — (l-B2)^ 
m 


cos  cp 


3G 

30 


cos 


0 sin  cp  3G  \ 
sin  0 3c p / 


- 1 


- n «V/2 


9G  + 
30 


2.51  p 1141  V2G 
6 


(l-02)3/2  6(r-r  *)  6(0-0')  6(0-0*)  6(c^cp*)  6(t-Q  . 

8 sin  0 V 

We  choose  the  initial  condition  G = 0 for  t > t^. 

Multiply  Eq.  (3.2)  by  0G  sin  0/ (1-0  ) 7 ; Eq . (3.5)  by  0W  sin  0/ 
2 3/2 

(1-0  ) ; add;  and  integrate  over  all  values  of  r,0,0,cp,t.  The  sev- 

eral terms  of  the  integration  on  the  left  lead,  after  integration  by 
parts,  to  the  following  boundary  terms: 


t=oo 

r=°°  1 

0=1 

GW 

, GW 

, GW 

, sin  0 GW 

t=-°° 

II 

o 

0=0 

Cp®=  2 7T 

sin  efe 

G3w\ 

0— 7T 

( W3G 

_ G3w\ 

cp2ir 

* 

cp=0 

\30 

~ 30  ) 

0=0 

\3cp 

‘ 3cp  / 

cp=0 

(3.6) 
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The  first  term  vanishes  at  t = 00  from  the  condition  on  G,  and  at  t = 
from  the  condition  on  W.  If  t is  finite,  W vanishes  for  r -*■  «>.  We  take 
G(r  = 0)  = 0.  Since  the  function  K(cos  0),  Eq.  (2.34),  has  an  upper 
limit,  there  will  be  no  electrons  at  6 = 1,  so  the  third  boundary  term 
vanishes  at  3 = 1.  To  make  the  lower  limit  contribution  vanish,  we 
take  G = 0 for  6 < 8*.  We  assume  G and  W are  bounded  functions  of  8, 
which  makes  the  0 boundary  terms  vanish,  and  that  they  are  periodic 
functions  of  the  velocity  azimuthal  angle  cp  with  period  2tt,  as  is  re- 
quired by  the  physics.  As  a result  the  cp  boundary  terms  vanish  and  the 
entire  left  side  integration  is  zero.  On  the  right  side,  the  second 
term  can  be  evaluated  completely  because  of  the  delta  functions,  and 
is  simply  -W(r  ' , 8^ , 0 ” , cp  ^ , t ') . Thus,  W can  be  expressed  in  terms  of  the 
Green's  function  G by: 

w(r',e\e',«p',o  = 


/“  /*°°  /*1  /*it  /*2tt 

dt  J dr  J d8  J J sin  0 d0dcp  f (t  - r/c)  K(cos  0)  6 (fi  - h(cos  8)) 


o(r,B,8,cp,t,r  ,8  ,8  ,cp*, 
2tt 


(3.7) 


The  time  may  be  eliminated  from  Eq.  (3.5).  It  has  the  following 
solution,  which  meets  the  initial  and  boundary  conditions: 


When  this  form  is  substituted  into  Eq.  (3.5),  the  terms  in  6'  cancel, 
and  the  terms  in  <5  factor  out  on  both  sides.  The  integral  in  Eq.  (3.8) 
can  be  evaluated,  and  the  argument  of  the  delta  function  becomes 
t-t ' + (y-Y*)/1.16  p.  Since  this  must  be  zero,  we  see  that  the  energy 
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of  each  electron  decreases  linearly  with  time,  exactly  as  it  must  for 
the  purely  coherent  energy  loss  process.  Substitute  Eq.  (3.8)  into 
Eq . (3.7)  and  evaluate  the  t integration,  and  we  have: 

W(r  ' ,8' ,Q' ,y' ,t')  = 


/*°°  /*!  /*27T 

/.  drX  de/o  L 


sin  9 d0d<p  f t' 


i + (y'-r) 

C 1.16  p 


K(cos  6)  6(8  - h(cos  8))  H(r  ,8,8,cp,r*,B'.9'.ccO 

(2tt  1.16  p) 


(3.9) 


The  differential  equation  satisfied  by  H will  be: 


9H 

98 


8 c cos  8 9H 
1.16  p(l-82)3/2  9r 


eB  8 / 9tl  cos  8 sin  cp  8h\ 

m 1.16  p(l-62)  V ? " ' Sin  0 9cp/ 


_ 2.164  2 6(r-Q  6(6-8")  6(8-0*)  6(<p-<p')  . 

Bd-B2)^  sin  6 


(3.10) 


We  shall  define  the  three  functions  q(8),  e(8)>  and  a (8)  as  the  coef- 
ficients in  this  equation: 


q (6)  = 


1.16  p(l-02) 3/2 


(3.11a) 


e(6) 


eB8 

m 1.16  p(l-82) 


(3.11b) 


a (6) 


2.164 

eu-e2)*5 


(3.11c) 


We  now  tackle  the  solution  of  Eq.  (3.10). 
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ELIMINATION  OF  THE  FIRST  DERIVATIVES 

The  procedure  for  simplification  of  Eq.  (3.10)  is  indicated  on 
p.  39  of  Ref.  11.  The  procedure  is  to  introduce  as  new  independent 
variables  a set  of  first  integrals  of  the  motion  in  the  absence  of 
scattering.  When  the  differential  equation  is  expressed  in  these  vari- 
ables, the  first  derivatives  will  disappear  from  Eq.  (3  10),  leaving 
only  3H/9B  and  second  derivatives  from  the  scattering  term.  The  new 
variables  are  essentially  the  integration  constants  of  the  motion  of 
an  electron  under  magnetic  deflection  with  energy  loss,  with  the  velocity 
magnitude  (3  regarded  as  the  independent  variable. 

In  the  absence  of  scattering,  the  y-direction  cosine  of  velocity 
will  remain  constant,  whereas  the  x-  and  z-direction  cosines  will  cor- 
respond to  a steady  rotation.  Expressing  the  direction  cosines  in  terms 
of  their  initial  values,  we  take  as  the  trial  representation  for  the  new 
variables 


sin  6 cos  cp  = sin  9 cos  cp  cos  £( 8)  - cos  0 sin  £(8)  (3.12a) 

o To  o 

sin  9 sin  cp  = sin  0^  sin  cpQ  (3.12b) 

cos  0 = sin  0 cos  cp  sin  £(B)  + cos  0 cos  £(B)  (3.12c) 

o T o o 

r = r - c sin  0 cos  cp  S(B)  + cos  0 C(8)  (3.12d) 

o I o To  o J 


This  expresses  the  old  variables  9,cp,r  in  terms  of  the  new  variables 

6 , cp  , r . The  functions  £(B),  C(B),  S(B)  are  to  be  determined  to  elim- 
o o o 

inate  the  first  derivative  terms  in  Eq.  (3.10).  Since  the  initial  value 
of  the  magnitude  of  the  velocity  is  S',  we  must  choose  £(8*)  = C(8  ) 

“ S (8  ) ■ 0.  The  inversion  of  Eq.  (3.12)  to  express  the  new  variables 
in  terms  of  the  old  yields: 

sin  9 cos  cp  • sin  0 cos  cp  cos  £(B)  + cos  0 sin  £(B)  (3.13a) 

o To 

sin  9 sin  cp 
o To 


sin  0 sin  cp 


(3.13b) 
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cos  0 = - sin  0 cos  cp  sin  5(8)  + cos  0 cos  £(0)  (3.13c) 
o T 

rQ  = r + c sin  0 cos  cp  S(B)  + cos  0 C(B)^  (3.13d) 
C(B)  = C(B)  cos  5(8)  + S (6)  sin  5(8)  (3 . 13e) 
S(B)  = C(B)  sin  5(0)  - S(0)  cos  5(B)  (3.13f) 


The  reduced  Green's  function  H will  be  denoted  by  I when  expressed 
in  terms  of  the  new  variables.  The  derivatives  of  H may  be  expanded  by 
the  chain  rule  in  terms  of  derivatives  of  I and  derivatives  of  the  new 
variables  with  respect  to  the  old.  The  latter  set  is  found  by  differen- 
tiating Eq.  (3.13).  The  manipulation  is  tedious  and  will  not  be  pre- 
sented. To  save  writing,  drop  the  arguments  from  5,C,S,  and  let  d/dB 
be  denoted  by  a prime. 

When  the  expressions  for  the  derivatives  of  r ,0  ,cp  are  substituted 

oo^o 

into  the  chain  rule  relations  to  obtain  the  derivatives  of  H,  and  the 
resulting  forms  are  substituted  into  Eq . (3.10),  an  enormous  amount  of 
algebraic  simplification  takes  place.  The  following  formulas  are  ob- 
tained for  the  coefficients  of  the  first  derivative  terms: 


sin  0 cos  cp ( S ^ + 2aS  - eC)  + cos  0(C^  + 2aC  + eS  - q) 


('O' 


(5*-e)  cos  cp  - a cos  0 
^ lo o 

sin  0 


(3.14a) 


(3.14b) 


ft)- 


(5'-e)  cos  0q  sin  cpo 
sin  0 


(3.14c) 


The  coefficients  of  sin  0 cos  cp  in  Eq.  (3.14a)  may  be  set  equal 
to  zero,  which  gives  a pair  of  differential  equations  for  C and  S which 
may  be  solved  easily.  Setting  5'  = e eliminates  Eq.  (3.14c)  and  the 
first  term  of  Eq.  (3.14b).  The  second  term  is  part  of  the  transformed 
Laplacian. 
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The  solutions  for  the  functions  C,C,S  satisfying  the  proper  initial 
conditions  are: 


-/  dBi  e(Bi’ 


(3.15a) 


•L  dBi  ’(B1 


-2  | «2a(B2) 

)e  1 cos  I d03  g(83) 

*'6l 


(3.15b) 
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(3.15c) 


= /r  d3l  q(6l 

"X-  d6i  q(6i 


-2  I d62  a(82)  f B, 

)e  1 cos  I 

Jr' 


cos  I dB,  e(BO  (3.15d) 

•V  J J 


-*/’ 
. *'B1 

)e  1 


d82  a(82) 


f 1 

sin  J dB3  e(83)  (3 . 15e) 


It  is  clear  from  these  equations  that  t,  represents  the  angle  through 
which  the  velocity  of  the  electron  would  turn  in  the  magnetic  field  while 
changing  velocity  magnitude  from  Q'  to  8 in  the  absence  of  scattering. 

The  functions  C and  S (C  and  S)  describe  the  forward  and  sideward  motion 
of  the  electron  during  this  turn,  and  the  exponential  is  associated  with 
the  effects  of  scattering. 

The  coefficients  of  the  second  derivatives  must  now  be  computed. 
These  require  straightforward  but  lengthy  manipulations.  The  result  is 
the  transformed  differential  equation: 


»T 
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6(rQ-r')  6(B-6')  6(0q-0')  6(cpo-cp') 
sin  0 


(3.16a) 


1 = 0 6 < B' 


(3.16b) 


APPROXIMATE  SOLUTION  OF  THE  DIFFERENTIAL  EQUATION 

Define  the  differential  operators  on  the  second  and  third  lines 

of  Eq.  (3.16a)  as  M and  Q respectively.  We  see  that  the  first  line 

differential  operator  includes  only  derivatives  with  respect  to  the 

initial  velocity  angles  0 ,cp  . The  second  line  differential  operator 

o To 

M involves  mixed  derivatives  with  respect  to  initial  position  r^  and 
initial  velocity  angles,  and  the  third  line  differential  operator  Q 
contains  only  position  derivatives. 

We  would  expect  that  the  primary  effect  of  scattering  for  small 
changes  in  velocity  magnitude  (early  times)  would  be  to  change  the 
velocity  direction.  As  the  velocity  change  increases,  the  effect  on 
position  becomes  significant.  We  would  therefore  anticipate  that  an 
expansion  in  terms  of  the  "ratio"  of  the  operators  M and  Q to  the  first 
line  operator  would  in  some  sense  be  rapidly  convergent  for  small  ve- 
locity changes,  and  hope  that  this  convergence  persists  for  the  signifi- 
cant duration  of  the  current.  We  shall  perform  this  expansion  to  third 
order  and  obtain  after  the  fact  that  the  convergence  has  remained 
satisfactory. 


We  see  from  Eqs.  (3.15d,e)  that,  for  small  changes  in  velocity 
from  the  initial  velocity  6",  C is  proportional  to  (8-6")  and  S to 

„ 2 

(8-6  ) • The  operator  M is  thus  a first-order  and  Q a second-order 
operator,  both  in  the  velocity  change  and  in  the  order  of  the  deriva- 
tive with  respect  to  r^.  Let  us  introduce  the  operator: 


L = 


ll  - 


a 

— J + cot 

30 

o 


9 

o 


3 . 2 

30~  + csc 

o 


(3.17) 


There  should  be  no  confusion  between  the  differential  operator  L and 
the  ionization  lag  function  L of  Eq.  (2.44b).  We  now  make  the  symbolic 
expansion: 

1=1  + I,  + I.  + I„  (3.18) 

o 1 2 3 

substitute  into  Eq.  (3.16a),  and  match  terms  of  corresponding  orders. 
There  results: 


6 0—0  6(8-8")  6(0o-8")  6 (cpQ— cp " ) 


LI  = ° 

u 

U 

(3.19a) 

o 

sin  0 

LI,  = MI 

(3.19b) 

1 o 

LI2  = MI1  + 

(3.19c) 

LI3  * MI2  + QI^ 

(3. 19d) 

The  expansion  can  be  continued  indefinitely,  but  the  calculation  of  the 
higher  terms  becomes  so  complex  that  we  have  not  attempted  to  proceed 
beyond  third  order,  and  have  in  some  parts  of  the  calculation  stopped 
at  second  order. 

From  the  structure  of  Eqs.  (3.19),  we  see  that  we  wish  to  solve 
the  general  equation: 
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Lu  = It  - a(6) 


»2  cos  0 „ „2 

d + iIirra^  + -TT7-7f  = R(6,eo,cpo)  (3.20) 

do  00  sm  0 dtp 

o o Yo_l 


Since  the  operator  L does  not  involve  r^,  the  dependence  of  the  right 
side  on  r^  need  not  be  indicated,  and  u is  to  be  regular  in  ®o»90* 

This  equation  is  obviously  suitable  for  an  expansion  in  Legendre  func- 
tions in  0o,cpo>  Thus,  we  assume  the  following  solution  form,  in  which 
the  orthogonality  properties  of  the  Legendre  functions  yield  the  co- 
efficient representation: 


'‘i 

u =]C  Sli  (i+1"2)  Pi(cos  eo)  [Di(6)  cos  jcPo  + EJ(6)  sin  Wo] 

• _ A • A 


i=0  j=0 


D-] 

1 /•  7T  /•2tt 


1 /•  TT  rm 

= | I sin  0.^  dO^dcp^  (cos  0 ) utBjO^cp..  ) 

1 Jo  Jo 


(3.21a) 

cos  jcp1 

(3.21b) 

sin  jcp1 

(3.21c) 

Here  we  have  written  = 1 if  j = 0,  2 if  j / 0,  and  have,  in 
anticipation  of  later  results,  changed  the  letter  for  the  dummy  inte- 
gration variable.  We  assume  a corresponding  expansion  for  R,  with  co- 
efficients Substitute  these  expansions  into  Eq.  (3.20),  multiply 

by  the  factors  of  Eqs . (3.21b,c),  and  integrate.  The  differential  oper- 
ator in  0O»9O  in  Eq.  (3.20)  reduces  to  -i(i+l)  when  applied  to  the 
appropriate  term  in  the  expansion.  The  orthogonality  of  the  Legendre 
functions  then  completely  separates  the  terms  of  the  expansion,  leaving 
the  differential  equations 


|g  dJ  + i(i+l)  a(B)  d|  = t| 


e]  + i(i+l)  a (6)  e|  = u] 


(3.22a) 


(3.22b) 
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Ass  inning  that  everything  vanishes  at  8 = 3',  the  solution  of  these 
equations  is: 


D1  -/»  deiTl  <6i> 


-i(i+l)  I d62a(62) 
e Jh 


(3.23) 


and  correspondingly  for  E'j.  Put  these  solutions  into  the  expansion 
(3.21a),  obtaining 


u(e’0o’cpo) 


e.  „ „ d82*<62> 

2tt  U+ 2)  (i+j)!  e 1 


i=0  j-0 


(cos  6g)  P^ (cos  0i)  cos  j(<po~cp1)  R(61,01,cp1) 


(3.24) 


Using  the  addition  theorem  for  Legendre  functions,  this  becomes  symboli- 
cally: 


u(8,eo,cpo)  = 

rB  f TT  /*2TT 

I d^  I I sin  0^^  dO^d^  R(B1,01,cp1)  X(8,0Q,cpo,61,01,(p1)  (3.25a) 

J&'  Vo  Jo 


X(6, 0Q  ,cpo 


» -i(i+l)  I dB2a(82) 

^ (i+4)e  ®1  P^cos  0q  cos  0^  + sin  0q  sin  0^  cos  (tp^-cp^)) 

i-0 

(3.25b) 
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The  function  X is  the  Green's  function  for  the  differential  equation 

(3.20).  We  shall  write  the  arguments  of  X symbolically  as  (£  ,£,), 

o 1 

where  the  letter  stands  for  the  triad  and  will  write  dto^ 

for  sin  0^  d0^dcp^.  In  all  these  equations,  &'  stands  for  a value  of 
3 immediately  below  B' . We  now  apply  this  general  solution  to  the 
sequence  of  Eqs.  (3.19a-d).  We  see  immediately  that: 

ID  = 6(ro-r')  X(£o.O  0.26) 


The  differential  operators  M,Q  may  be  split  into  products  of  9/3rQ 
and  operators  in  8,0Q,cpo,  which  we  denote  M^,Q^.  The  9/3rQ  will  act 
to  produce  derivatives  of  6(r  -r'),  whereas  the  other  operators  will 
appear  in  integrals  of  the  type  in  Eq.  (3.25a).  We  shall  write  down 
the  results  of  the  successive  solution  of  Eqs.  (3.19a-d)  without  the 
intermediate  details: 


1 = 6(r  -r  *)  X(f,  .O  - 2c  6 
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(3.27) 
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This  representation  for  I holds  for  g > g'.  For  g < g',  1=0.  At 
g = g^,  I reduces  to  the  first  term.  If  we  let  g approach  from 
above,  X tends  to  6(0o~0')  6 /sin  0,  so  I satisfies  the  proper 
discontinuity  condition. 

The  reduced  Green's  function  H is  found  from  1 by  re-expressing 
the  variables  r0>90»90  in  terms  of  r ,0,<p,  using  Eq.  (3.13).  Then  H is 
to  be  substituted  in  Eq.  (3.9)  to  obtain  W,  the  response  to  the  Compton 
production  function,  as  represented  in  the  variables  This 

representation  is  then  to  be  inserted  into  Eqs.  (3.3)  and  (3.4)  and, 
after  a change  of  lettering,  integrated  over  E,' . The  result  will  be 
the  current  and  the  ionization  production  rate.  We  shall  carry  out 
these  procedures  che  next  section. 
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IV.  CALCULATION  OF  THE  CURRENT  AND  IONIZATION 


The  highest  order  term  in  the  expansion  (3.27)  involves  integra- 
tion over  nine  variables.  There  are  four  integrations  in  Eq.  (3.9) 
and  three  in  Eq.  (3.3)  or  (3.4).  Thus,  we  have  to  deal  with  16-fold 
integrals.  We  shall  now  perform  analytically  as  many  as  we  have  found 
possible.  We  shall  interchange  orders  freely  but  carefully  and  will 
make  great  use  of  the  orthogonality  properties  of  the  Legendre  functions. 

The  first  integration  is  the  outermost,  the  integration  over  &',cp, 
in  Eqs.  (3.3)  and  (3.4).  We  observe  that  in  Eq.  (3.27)  these  variables 
occur  only  in  the  extreme  right  factors  of  each  term  in  the  expansion. 

In  X (£,£'),  the  angular  variables  appear  only  as  arguments  of  Legendre 
functions.  We  thus  have  to  evaluate  the  integrals: 


^cos  0^  cos  Q'  + sin  0^  sin  0*  cos  (cp^-’P^J  (4.1a) 


/ 


dw'  sin  0'  cos  cp 


(4.1b) 


where  U^cp^  refer  to  the  first  arguments  of  X in  the  appropriate  term 
of  Eq.  (3.27).  Rotate  the  coordinate  system  so  that  the  pole  is  in  the 
direction  0^,tp^,  and  Eq.  (4.1)  becomes: 


(cos  0”) 


J dw'  PJ 

/*"'[  sin  0^  cos  cp^  cos  0'  + cos  0^  cos  cp^  sin  0'  cos  cp' 
- sin  sin  0"  sin  P^(cos  0") 


(4.2a) 


(4.2b) 


The  first  integral  is  zero  unless  i = 0,  in  which  case  it  equals  47T . 
The  second  and  third  terms  in  the  bracket  of  Eq.  (4.2b)  vanish  when 
integrated  over  cp'.  The  first  integrates  to  zero  unless  i = 1,  when 
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it  becomes  4tt  sin  0^  cos  9^/3.  Using  these  results  in  the  series  repre- 
sentation then  yields: 


f 


du>'  x(Sk.O  = 1 


/ 


do)'  sin  0"  cos  cp'  = sin  0^  cos  cp^  e 


-2  | dB~a(B~) 

6' 


(4.3a) 


(4.3b) 


Since  is  a differential  operator,  when  it  is  applied  to  a con- 
stant the  result  is  zero.  Hence,  there  is  no  contribution  to  the 
ionization  density  from  the  second,  third,  fifth,  and  seventh  terms  of 
the  expansion  (3.27).  We  shall  make  a further  simplification  by  drop- 
ping the  third-order  terms  in  the  current.  This  will  be  justified  by 
later  numerical  work. 

Next  are  the  argular  integrations  in  Eq.  (3.27).  We  shall  demon- 
strate the  procedure  by  calculating  the  first-order  term  in  the  current 
and  then  writing  out  the  result  for  the  complete  expression.  In  this 
term,  the  integral  to  be  evaluated  is: 


/ 


du^  X(?o,C1)  M1(^1)  sin  0^^  cos 


(4.4) 


The  differential  operator  is  given  in  Eq.  (3.16a),  with  the  variables 
to  be  given  the  subscript  1.  There  results: 


M^(f^)  S*n  ^1  COS  ^1  = C1  s*n  ®i  cos  ®i  cos  9^ 


- S^l  - sin2  0^  cos2  9^) 


(4.5) 


Insert  this  into  the  integral  in  Eq.  (4.4)  corresponding  to  the  itTi 
term  in  the  series  for  X,  and  rotate  the  coordinate  system  so  that  the 
pole  is  in  the  0O»9O  direction.  We  then  have,  after  simplification: 
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P^(cos  0^)  SL  + (sin  0q  cos  cpQ  cos  0^  + cos  0q  cos  cpo  sin  0^  cos  <p^ 

- sin  m sin  0 sin  cp.  ) < C, (cos  6 cos  0. 

o i I ( 1 o 1 

- sin  0 sin  0,  cos  cp. ) + S,  (sin  0 cos  c p cos  0, 

o 1 T1  1 o ~o  1 

+ cos  0^  cos  cpQ  sin  0^  cos  cp^  - sin  cp^  sin  0^  sin  cp^)jj 


(4.6) 


Integrating  over  cp^  yields: 


2tt  f sin  0^  d0^P^(cos  0^) 

'o 


cos  0 sin  0 cos  cp  [ 4 cos2  0,  - ^ ) 
1 o o Ho\  2 12/ 


+ j - 1 + sin2  0^ 


2 2 
+ sin  0 cos 
o 


cpo(|  cos2 


(4.7) 


From  the  orthogonality  properties  of  the  Legendre  functions,  the  inte- 
gral is  zero  unless  i = 0 or  2.  Evaluating  these  cases  and  restoring 
the  series  representation  gives  the  result: 


f^l  x(C0.?i)  M^f^)  sin  0^^ 


cos 


Ti 


-6  / d£La(6~) 

8 ^ L 

C-  cos  0 sin  0 cos  cp  e * P1 
1 o o To 


f 

Je 


6 


- s. 


l + /l  _ 

3 ^3 


*•) 


-6  / dg .3(6  ) 

. 2 - 2 \ J Q,  2 

sin  0 cos  <p  I e 1 
o 


6 


(4.8) 
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In  this  expression  we  have: 


and  a corresponding  expression  for  S^.  All  the  angular  integrals  of 
Eq.  (3.27)  can  be  evaluated  in  exactly  the  same  manner,  but  the  manipu- 
lations become  very  tedious  for  the  higher  order  terms.  We  write  out 
the  result  to  third  order  in  the  integrand  for  the  ionization  and  to 
second  order  in  the  current: 


/' 


I dm  * * 6(r  -r  *)  4-  c 


2 6 lro‘r’)  fr  d6la(6l>  (C1  + Sl)(j  + 3 


/•B 

-6  / d6,a(6,> 

e ‘ 


'(■ 


. cos  0 4-  S,  sin  0 

1 o 1 


cos  i 
o o I 


-6 y B dB2a(82) 


2 71). 

e l 


/ B,  -b  I dBjatBj) 

dB,a(6.)  / dB,a(B,)e  ’,82 
8'  1 JB'  22 


(C1C2  + S1S2)(C2  cos  ®o  + S2  sln  6o  COB  'P0)(2C 


r8  /-B 

-12  / d6.a(6.)  -2  / dfl  . 


^(8,) 


cos  0 4-  S,  sin  0 cos  ) | e 1 

o 1 o n • 1 


+ (C2  + S0(C1  “* 

- 5 [c.  cos  0 4-  S.  sin  0 cos  ] ( C_  cos  0 4-  S~  sin  0 cos  cp  1 

V 1 ° 1 ° °)  \ 2 ° 2 ° °/ 


-12  dBja(B))  -2 /;  dB3a(B3 

/, 


-12  I dB.alB,) 
2 2B, 

e 1 


(4.10) 
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The  next  step  is  to  replace  the  variables  r0>®0*'P0  by  r,0,cp  using 
Eq.  (3.13).  The  integrations  which  remain  to  evaluate  the  current  and 
the  ionization  density  are  over  the  variables  r,8,9,<p,B^,B^,B2.  The 
variable  r appears  in  the  argument  of  a delta  function,  r-r'  + c 
(C  cos  0 - S sin  0 cos  <p) , its  derivatives,  and  f.  The  derivative  with 
respect  to  r of  the  delta  function  goes  over  into  3/9(ct')f  after  the  r 
integration,  where  the  argument  of  f is  t'-r'/c  - (y-y')/1.16  p + 

C cos  0 - S sin  0 cos  cp.  Since  t'-r'/c  is  the  retarded  time  at  the 
electron  location,  we  have  now  shown  that,  to  the  approximation  of 
neglecting  the  variation  of  p and  g over  the  electron  motion,  the  cur- 
rent and  ionization  are  functions  of  retarded  time  only,  as  is  of  course 
to  be  expected. 

The  variable  6 appears  as  the  argument  of  a delta  function  8 - 
h(cos  0)  and  in  integration  limits.  Since  8 must  exceed  the  &' 
integration  must  be  brought  inside  the  8 integration,  and  then  the 
limits  on  become  0,h.  Referring  back  to  Eqs.  (3.15b,c),  when  8 is 
replaced  by  h,  we  shall  regard  C,S  as  functions  of  the  lower  limit  &' . 

We  then  define: 


= c(B')  - C(B1)  = 


-2 

d$2q(&2)e 


d63a(83) 


d84e(84) 


(4.12) 


and  a corresponding  equation  for  S^.  The  function  C(&')  differs  from 
Eq.  (4.12)  only  in  that  the  outermost  upper  limit  is  h.  We  have  the 
relations: 


C 


1 


S 


1 


d82a(B2) 

(Cx  cos  X,  + S sin  t) 

d82a(B2) 

(C^  sin  C - S^  cos  ?) 


(4.13a) 


(4.13b) 
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C-  cos  0 + S.  sin  0 cos  c p = e P1 

1 o 1 o Yo 


2Je  d*2a<62> 


(C^  cos  0 - sin  0 cos  (p) 


(4.13c) 


2 2 

= e 


d82a(62) 


(c2  + s2) 


(4.13d) 


2 JR  d63a(B3>  + 2 JR  d*3a(B3> 

C1C2  + S1S2  = e *2  ^—1—2  + -1-2J  (4 . 13e) 


The  point  of  introducing  these  functions  is  that  C(B  ) and  S(B^),  func- 
tions of  only  one  variable,  are  easy  to  compute.  The  underscored 
functions  are  then  obtained  by  subtraction,  or  the  polynomials  in 
C(B^)  - C(0j)  may  be  expanded  and  the  integrations  thereby  separated 
into  products  of  integrals.  This  procedure  greatly  shortens  the  com- 
puting time. 

These  discussions  show  how  to  integrate  over  the  variables  r,0. 

We  shall  now  assume  that  f(t  - r/c) , the  gamma-ray  time  history,  is  a 
delta  function  whose  argument  is: 


r 3d 8 ],  n J 6 

t - r/c  -I  r 2~3/2  ) 1 " 6e 

6 1.16  p (i-8  ; i 1 } 


-2  f »<«!)«! 


|cos  0 cos  J d82e(82)  - sin  0 cos  cp  sin  dB2e(82)j 


(4.14) 


The  derivatives  of  the  delta  function  with  respect  to  t-t^  which  arise 

from  the  r derivatives  in  I may  be  taken  outside  the  integrations, 
o 
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Thus,  Eq.  (4.14)  indicates  how  we  may  integrate  over  8^.  Choose  values 
of  t'-r'/c,  9.  and  <p.  Find  the  value  of  S'  such  that  Eq.  (4.14)  is 
zero.  Set  8^  equal  to  that  value  in  Eqs.  (4.10)  and  (4.11).  Divide 
by  the  derivative  of  Eq.  (4.14)  with  respect  to  S',  which  is  of  course 
the  integrand  of  Eq.  (4.14).  The  result  will  be  the  integral  over  S' • 

The  structure  of  Eq.  (4.14)  is  such  that  the  integral  increases 
monotonically  as  S'  decreases  from  h.  At  8"  = h,  expression  (4.14)  is 
t'-r'/c  and  is  positive.  When  S'  = 0,  the  integral  has  a maximum 
value  which  is  a function  of  0 and  cp.  For  sufficiently  small  values  of 
t'-r'/c , there  will  be  a zero  for  Eq.  (4.14)  as  a function  of  8 for 
all  values  of  0 and  cp.  As  t'-r'/c  increases,  it  reaches  a value  which 
exceeds  the  maximum  value  of  the  integral  for  some  values  of  0,cp,  which 
then  do  not  contribute  to  the  integration  over  S' . Those  values  of  0,cp 
correspond  to  electrons  which  have  slowed  to  zero  velocity,  so  of  course 
they  should  not  contribute.  Eventually,  t'-r'/c  becomes  so  large  that 
it  exceeds  the  integral  for  all  0,cp  and  there  is  no  further  contribution 
to  the  current  or  ionization. 

We  have  reached  this  point  in  the  solution  of  the  transport  equa- 
tion and  the  calculation  of  the  current  and  ionization  without  utilizing 
the  specific  forms  of  the  functions  q,  e,  and  a.  We  shall  now  employ 
Eq.  (3.11)  to  simplify  the  expressions.  We  first  evaluate  £ and  the 
integral  of  a.  With  Eqs.  (3.11b)  and  (3.15a),  and  8 set  equal  to  h, 
we  have: 


fh  /*h  B.dB, 

■ / dB  e(B  ) = ~ Z f ~ — Y~ 

J S'  1 m 1.16  p J S'  (1-Bp 


(4.15a) 


eB  1 


1.16  p 


— log 


VS 


(4.15b) 


1 i Y° 

k log  r 


(4.15c) 
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where  the  constant  K is  equal  to: 


K = 6.595  £ 

D 


(4.16) 


and  Yq  is  the  initial  electron  total  energy,  including  the  rest  energy 
in  units  of  mc^. 

We  can  invert  Eq.  (4.15c)  to  express  B^  and  y'  in  terms  of  , 
yielding: 


(4.17a) 


r - i - 


2K5  * \ 4 


(4.17b) 


Thus,  as  the  electron  turns  in  the  magnetic  field  through  an  angle  X, 
in  the  absence  of  scattering,  the  energy  decreases  exponentially  with  £. 
Since  the  smallest  possible  value  of  B is  zero,  we  see  from  Eq . (4.15c) 
that  the  maximum  angle  of  turn  of  the  electron  is  directly  proportional 
to  the  magnetic  field,  inversely  proportional  to  the  density,  and  pro- 
portional to  the  logarithm  of  the  initial  energy.  For  electrons  with 
the  energy  distribution  of  Fig.  7,  the  logarithm  has  a maximum  value 
1.306  at  0 = 0 and  is  down  to  0.7  at  40  deg.  We  shall  consider 
log  Yq  = 1 as  representative,  corresponding  to  electrons  slightly  below 
the  mean  energy,  0.8  MeV.  With  a magnetic  field  of  0.6  gauss,  the  maxi- 
mum angle  of  turn  will  be  72  deg  at  an  altitude  of  20  km,  347  deg  at 
30  km,  and  1568  deg  (4.36  revolutions)  at  40  km.  This  cutoff  of  the 
angle  of  turn  is  caused  entirely  by  the  energy  loss  process. 

The  integral  of  a can  be  evaluated  using  Eq.  (3.11c).  We  shall 
arrange  the  limits  in  the  exponents  of  Eqs.  (4.10)  and  (4.11)  so  that 
the  upper  limit  is  always  h.  After  simplification,  we  obtain: 


r -“W  ■ 1-082  ^ 

J p L o 


(4.18a) 
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(Y  - 1)  (Yo 
(Y  + i)  (y0 


+ 1) 

~=~rj 


2.164 


(4.18b) 


We  shall  introduce  ( as  a new  independent  variable.  We  have,  for 
the  function  C (see  Eq.  (4.12)),  the  expressions: 


c(«'>  "Sb 


£'«(<: 


-2K£ 


y eK^-ly  +1 
o __  1 o 

Y e'K?  + 1 Yo  " 1 
o 


2.164 


cos 


c 

(4.19) 


For  S,  cos  X,  is  replaced  by  sin  X,.  This  representation  depends  directly 
on  X,' , related  to  8"  through  Eq.  (4.17b),  and  indirectly  on  9,  which 
appears  in  y . Equation  (4.19)  is  in  precise  form  for  step-by-step 
numerical  integration. 

If  there  were  no  scattering,  we  would  set  a = 0.  The  functions 
C,S  then  represent  the  forward  and  sideward  components  of  position  of 
an  electron  with  initial  energy  Yc  which  has  turned  through  an  angle 
X,'  from  its  initial  direction.  The  first  effect  of  scattering  is  to 
insert  the  factor  (4.18b)  into  C and  S.  This  corresponds  to  a change 
in  the  velocity  of  a mean  position.  If  we  were  to  keep  only  the  first 
terms  in  the  expansions  (4.10),  the  motions  of  the  aggregate  of  elec- 
trons would  be  fully  described  by  a moving  mean  position. 

The  scattering  expression  (4.18a)  corresponds  closely  to  Longmire's 
obliquity  theory,  which  sets  the  right  side  of  (4.18a)  equal  to  (n— 1) /2 
(see  Eq.  (22)  of  Ref.  8,  with  a slight  change  in  constants).  Thus,  in 
some  sense,  the  obliquity  theory  is  a zero th  approximation  to  the  theory 
we  have  developed  in  this  report,  with  the  electron  velocity  distribu- 
tion Green's  function  reduced  to  the  first  term  of  I.  We  have  not 
attempted  to  work  out  the  full  details  of  the  comparison. 

We  shall  now  write  out  the  full  expressions  for  current  and  ioniza- 
tion density  as  integrals  over  Let  T = t'-r'/c,  the  re- 

tarded time.  We  define  the  integral  of  Eq.  (4.14),  the  equivalent  of 


■’ — — 


' 
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( h 

the  arrival  time  of  the  electron,  as  ^(c",6,<p) . Then  the  condition  that* 
Eq.  (4.14)  vanishes  determines  a value  such  that: 


m f*.'  /v  , Y + 1 \ 2.164 

T ■ T[C'(T,0.9),0.,1  Yd? 


(cos  0 cos  C _ sin  0 cos  c p sin  £)  (4.20a) 


T^C',0)  + T2(?',0)  cos  <p  (4.20b) 


The  value  of  X,'  to  be  used  in  the  current  and  ionization  is  to  be  found 
by  inverting  Eq.  (4.20).  In  practice,  this  must  be  done  by  interpola- 
tion. The  unspecified  limits  on  0 and  cp  are  for  all  values  such  that 
Eq.  (4.20)  has  a solution.  We  define  U(£)  as  the  function  of  Eq.  (4.18b), 
which  we  can  call  the  scattering  reduction  factor,  and  will  write  U'  for 
U(0 , for  U(£^),  and  U f°r  • A function  D is  defined  by: 


D(C>j9.cp)  = 1 - 8^U^(cos  0 cos  X,'  - sin  0 cos  cp  sin  t,')  (4.21) 

D is  the  ratio  between  retarded  time  and  local  time  of  the  aggregate 
of  electrons.  We  also  define  a function  A(£)  = a($)  d6/dC  as: 


A(0 


2.164  Ky 

Y2  - 1 


(4.22) 


We  introduce  the  shorthand: 

w = cos  0 sin  X,'  + sin  0 cos  cp  cos  X, ' (4.23a) 

Z.  = C.  cos  0 - S,  sin  0 cos  cp  (4.23b) 

l — i — i T 


We  now  have  the  final  result  for  the  current  and  ionization  rate: 


(4.24) 


J 


(4.25) 
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The  time  derivatives  in  these  equations  should  be  placed  outside 
the  integral  sign,  since  the  limits  are  functions  of  T.  We  have  written 
them  inside  to  avoid  having  to  repeat  the  first  line  of  Eqs.  (4.24)  and 
(4.25)  for  each  term,  but  in  calculation  the  time  derivatives  must  be 
treated  last. 

This  is  as  far  as  we  have  been  able  to  proceed  analytically.  We 
have  written  a computer  program  to  calculate  the  current  and  ioniza- 
tion from  Eqs.  (4.24)  and  (4.25).  This  program  is  purely  an  exercise 
in  numerical  integration,  since  all  the  differential  equations  have 
been  solved.  We  next  discuss  the  construction  of  the  program. 


! 


V.  DEVELOPMENT  OF  THE  COMPUTER  PROGRAM 


Although  the  programming  of  Eqs.  (A. 24)  and  (4.25)  is  fairly 
straightforward,  there  are  some  rather  tricky  numerical  points,  mostly 
associated  with  the  zero-velocity  cutoff.  The  program  takes  as  inputs 
the  magnetic  field,  gamma-ray  energy,  time-step,  and  angle-step.  Most 
calculations  have  been  done  with  a time-step  of  0.1  shake  and  an  angle- 
step  of  1 deg  in  the  turn-angle  C, . Double  precision  arithmetic  has 
been  employed  throughout. 

The  program  uses  different  numerical  integration  routines  for  the 
several  variables.  The  integration  over  0 is  a simple  summation,  since 
the  end-point  terms  of  the  trapezoidal  rule  vanish  because  of  the  factor 
sin  0.  We  have  used  an  integration  step  in  0 of  3 deg  for  0 below  30  deg 
and  6 deg  above.  This  is  a good  match  to  the  rate  of  variation  of 

K(cos  0)  sin  0 (see  Fig.  6).  The  innermost  integrals  over  £ and  use 

(17)  1 Z 

the  modified  trapezoidal  rule.  This  is  a four-point  rule,  which 

must  be  specially  modified  at  the  end  intervals  of  the  integration  range 
or  if  there  are  fewer  than  four  points  available  (as  can  occur  for  low- 
energy  electrons) . The  rule  reads 


f 


f(x)dx  = A*t9f(0)  + 19f  (1)  - 5f  (2)  + f (3)  1 

24 


(5.1a) 


r 

lAx 


(i+1)  Ax 


ffaOd,  ■ a-*.l.-_Ui-1>  13.ft.iU)..-  f (1+211  (5.lb, 


/, 


NAx 


(N-l)Ax 


f(x)dx  . WOir.l),  - ,5LLN_-_2L±  19  f (N-l)  + 9f(N)]  (5>lc) 


where  the  argument  lAx  of  f(lAx)  has  been  abbreviated  i and  similarly 
for  all  other  terms.  Equation  (5.1a)  represents  the  first  interval  of 
the  integration  range,  (5.1c)  represents  the  last  interval,  and  it  is 
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assumed  that  N is  at  least  3.  The  error  term  in  Eq.  (5.1b)  is 
5 (4)  (4) 

ll(Ax)  f (£)/720,  where  f'  ' (£)  denotes  the  fourth  derivative  of  f 
evaluated  at  some  point  in  the  interval.  The  functions  we  have  to  deal 
with  are  sufficiently  smooth  that  the  error  term  can  be  neglected 
throughout . 

The  integration  over  c p is  performed  by  gaussian  integration  pro- 
cedures. First,  we  reconsider  the  limits  of  integration.  Referring 

to  Eq.  (4.20),  we  define  a function  T (0,cp)  by: 

m 


V6’<P>  = IB 
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(5.2a) 


2.164 


Yo  = 


[(1+k) 2 + k2  cos2  0] 
[(1+k)2  - k2  cos2  0] 


(5.2b) 


This  is  the  time  for  an  electron  to  slow  to  zero  velocity  from  the 

initial  energy  Yq  appropriate  to  the  direction  0,cp.  It  is  clear  that 

T can  be  split  into  the  form: 
m r 


TJ0.'?)  = Tmi(0>  + Tm9(6)  cos  cp 


m2 


(5.3) 


This  lies  between  the  limits  T , ± T „ . Suppose  now  we  choose  a value 

ml  m2  r 

of  0,  and  let  T increase  from  zero.  As  long  as  T is  less  than  T , - T , 

ml  nw 

there  will  be  a solution  of  Eq.  (4.20),  T = T(£',0,cp)  for  all  values  of 

cp.  If  T is  between  T , - T . and  T , + T _,  there  will  be  a solution 
r ml  m2  ml  m2 

for  cos  cp  S (T  - anc*  not  otherwise.  If  T exceeds  T^  + 

the  integral  is  zero.  Since  the  integrand  depends  only  on  cos  cp,  the 
integral  from  0 to  2tt  is  twice  the  integral  from  0 to  TT.  We  thus  have 
the  integration  forms: 


/ 


f (cos  cp)dcp, 


L 


% 


f (cos  cp)  dcp 


(5.4a) 


whereas  the  second  integral  leads  to  a modified  version  of  Eq . (5.5): 


Here  n is  the  order  of  the  gaussian  approximation.  We  have  tested 
the  program  with  n = 5 and  n = 9.  The  results  were  indistinguishable 
to  three  significant  figures,  so  we  used  n = 5 to  reduce  the  running 
time . 

Other  difficulties  at  the  low-speed  range  are  caused  by  the  func- 
tion U.  As  can  be  seen  from  its  definition  (Eq.  (4.18b)),  near  the  cut- 
off turn-angle  £ = log  y /K,  U vanishes  like  (C  - inverse 

mo  m 

cubes  of  U appear  in  Eqs.  (4.24)  and  (4.25)  multiplied  by  combinations 
of  C's  and  S's  which  make  the  integrals  convergent.  However,  the  calcu- 
lation requires  division  of  small  quantities  by  small  quantities,  with 
attendant  loss  of  accuracy.  To  avoid  this  problem,  we  cut  off  the  £ 
integrations  one  step  before  the  last  interval,  and  then  extrapolate 
to  obtain  the  last  value.  Testing  has  shown  this  procedure  to  be  suf- 
ficiently accurate. 

With  these  mathematical  problems  clarified,  we  describe  the  general 
nature  of  the  program.  The  integration  rule  (Eq.  (5.1))  has  been 

Reference  17,  formula  25.4.38,  p.  889. 
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packaged  as  a subroutine  (Part  7)  which  integrates  many  functions  sim- 
ultaneously. The  atmospheric  density  p is  stored  for  altitudes  between 
18  and  60  km  in  2 km  steps.  Above  60  km,  the  current  per  primary  gamma 
is  essentially  independent  of  altitude,  and  the  ionization  is  negligible. 
The  program  can  take  any  altitude  for  which  the  density  has  been  stored 
as  an  input  initial  altitude. 

The  program  begins  at  the  initial  altitude  (18  km  if  not  otherwise 
specified)  and  the  first  value  of  0 (3  deg).  It  calculates  K(cos  9) 
and  h(cos  0)  = Using  Part  7,  the  functions  C,  S,  T^,  and  T^  (see 

Eqs.  (4.19)  and  (4.20))  are  found  and  stored  in  Part  2 for  all  £,  in 
1 deg  steps,  up  to  either  £ = log  or  £ = 120  deg,  whichever  occurs 

first.  The  value  120  deg  was  found  to  cover  the  full  time  range  of 
0 to  10  shakes  for  all  electrons.  Next,  the  functions  which  appear  as 
single  integrals  over  £ in  Eqs.  (4.24)  and  (4.25)  are  calculated  in 
Part  3.  The  expressions  C^,S^,  and  their  powers  which  involve  values 
of  C,S  at  two  values  of  Z,  are  expanded  using  Eq . (4.12),  the  separate 
single-value-of-C  functions  integrated  via  Part  7,  and  then  the  various 
functions  recombined  to  form  the  indicated  integrals.  This  procedure, 
which  permits  direct  step-by-step  integration,  is  much  faster  than  form- 
ing the  two-values-of-t,  functions,  which  requires  multiple  cycling.  (If 
there  are  N values  of  £,  and  the  relative  number  of  functions  after  ex- 
pansion is  m,  the  computing  time  is  reduced  by  a factor  N/2m.  Since  m 
is  about  2,  for  N = 120,  the  usual  value,  this  part  of  the  computation 
is  speeded  up  by  a factor  of  30.)  In  Part  4,  the  double  integrals  of 
Eqs.  (4.24)  and  (4.25)  are  calculated  by  the  same  procedure,  which  this 
time  uses  the  already  calculated  single  integrals,  come  of  which  appear 
as  the  inner  integrands  in  the  double  integrals.  At  this  point,  the 
factors  in  brackets  in  Eqs.  (4.24)  and  (4.25),  exclusive  of  the  time 
differentiations,  have  been  calculated  and  stored  for  all  £. 

In  Part  5,  the  time  T is  set  equal  to  the  first  time  step,  and  i 
(see  Eqs.  (5.5)  and  (5.6))  set  equal  to  1.  Equation  (4.20)  is  solved 
by  inverse  interpolation  to  find  the  value  of  Z'  to  be  inserted  into 
Eqs.  (4.24)  and  (4.25).  The  value  is  used  for  direct  interpolation  to 
find  the  values  of  the  terms  in  the  bracketed  expressions  of  Eqs.  (4.24) 
and  (4.25),  which  are  then  multiplied  by  the  exterior  factors  B^U^/D  or 
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1/D.  (The  factor  2tt  disappears — the  2 from  halving  the  c p integration 
range  and  it  from  its  use  in  Eq.  (5.5)  or  (5.6).)  The  value  of  i is 
iterated  from  1 to  5 and  the  expression  summed  using  Eq.  (5.5)  or  (5.6), 
whichever  is  appropriate  to  the  selected  values  of  0 and  T.  The  time 
is  stepped  up  to  either  11  shakes  (the  extra  shake  is  provided  for  lee- 
way in  later  differentiations),  or  to  the  last  value  such  that  there  is 

contribution  to  the  integral  (T  , (0)  + T „(0)),  whichever  occurs  first. 

ml  mz 

The  time  derivatives  are  then  calculated  using  five-point  Lagrange  dif- 
* 

ferentiation.  The  results  are  printed  and  stored.  There  are  many 
delicate  refinements  in  these  procedures,  mostly  associated  with  the 
transition  from  Eq.  (5.5)  to  Eq.  (5.6)  and  the  angle  and  time  cutoffs. 

The  angle  0 is  then  stepped  to  its  next  value  and  the  calculation 
returned  to  the  beginning  of  Part  2.  After  calculation  through  Part  5, 
the  results  are  summed  with  the  previous  expressions  to  form  the  inte- 
gration over  0.  The  successive  values  are  weighted  by  sin  0 K(cos  0) 
and  the  variable  value  of  the  0 increment.  This  process  is  continued 
to  the  last  value  of  0 (84  deg),  at  which  point  we  find  the  current 
and  ionization  rate  as  functions  of  time  at  the  indicated  altitude. 

The  ionization  rate  is  then  inserted  into  Eq.  (2.44)  and  integrated  in 
Part  6 to  obtain  the  total  ionization.  A special  integration  routine 
is  required  here  because  Eq.  (2.44)  cannot  be  integrated  stepwise  but 
must  be  cycled  in  t to  take  account  of  the  explicit  appearance  of  t in 
the  integrand.  Finally,  the  altitude  h is  incremented  and  the  entire 
procedure  repeated  until  the  highest  altitude  is  reached.  The  mean 
running  time  is  about  70  sec  per  altitude  step  on  the  IBM  370/158. 

This  completes  the  description  of  the  computer  program,  which  is 
included  as  an  appendix  to  this  report.  We  have  executed  the  program 
for  the  representative  values  B = 0.6  gauss,  E^  = 1.6  MeV,  and  shall 
present  and  discuss  the  results. 

Reference  17,  formula  25.3.6,  p.  883. 


VI.  RESULTS  AND  CONCLUSIONS 

We  shall  first  discuss  the  convergence  of  the  series  of  integrals 
in  Eqs.  (4.24)  and  (4.25).  In  Figs.  8a,  8b,  and  8c,  we  have  plotted 
the  ionization  per  primary  gamma  ray  versus  time  at  altitudes  of  20, 
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30,  and  40  km.  The  quantity  actually  plotted  is  10  Ng(x)/g(r),  and 

the  three  figures  have  different  vertical  scales.  The  terms  corre- 
sponding to  time  derivatives  of  various  order  in  Eq.  (4.25)  are  denoted 
by  subscripts,  thus  is  the  contribution  to  the  total  ionization  from 
the  second  derivative  in  Eq . (4.25).  The  curve  N is  the  sum  of  the 
three  contributions.  In  Figs.  9a,  9b,  and  9c,  we  have  corresponding 

g 

plots  for  10  ycJc(t)/2,  designated  as  J on  the  figures. 

First  consider  Fig.  8,  the  ionization.  It  is  apparent  that  the 
term  Nq  is  by  far  the  major  contributor.  Numerically,  at  20  km,  Nq  is 
at  least  84.5  percent  of  N,  the  minimum  occurring  at  T = 1.  The  second 
derivative  term,  N^,  makes  a maximum  contribution  of  9 percent  at  T = 1, 
and  provides  a maximum  of  7 percent  at  T = 0.5.  At  30  km,  Nq  is  at 
least  95.9  percent,  at  most  3 percent,  and  N ^ at  most  1.2  percent. 

The  time  integration  of  Eq.  (2.44a),  combined  with  the  time  differenti- 
ations of  Eq.  (4.25),  makes  the  contributions  of  N,,  and  N^  approach  zero 
at  large  T.  Figure  8 displays  very  clearly  the  rapid  convergence  of 
the  expansion  used  in  Eqs.  (3.19a-d)  as  far  as  ionization  is  concerned. 

The  convergence  of  the  expansion  for  the  current,  presented  in 
Fig.  9,  is  less  good,  but  still  satisfactory.  At  the  lowest  altitude, 

20  km,  the  early  peak  of  the  current  is  raised  and  sharpened  by  the 
higher  order  terms,  and  the  curve  decreases  more  rapidly  thereafter. 

The  same  effect,  only  less  pronounced,  occurs  at  the  higher  altitudes. 

At  20  km,  the  higher  order  terms  can  be  as  much  as  40  percent  of  the 
total  near  the  peak,  but  the  qualitative  behavior  is  not  changed  sig- 
nificantly from  Jq  and  there  is  no  physical  reason  to  suspect  major 
contributions  from  still  higher  order  terms.  In  any  case,  the  conver- 
gence of  the  series  is  somewhat  obscured  by  the  coupling  between  the 
time  derivatives  and  the  delta  function  gamma  source.  If  a smooth  source 
were  used,  the  series  would  converge  much  better. 


Contributions  to  the  ionization;  B=0.6  Ey=1.6MeV  H = 20  km 


Fig.  9b — Contributions  to  the  current;  B=0.6  E = 1 .6  MeV  H - 30  km 


-66- 


It  is  of  interest  to  see  what  effect  the  Klein-Nishina  distribu- 
tion of  initial  velocities  has  had  on  the  current.  In  Figs.  10a,  10b, 
and  10c,  we  have  plotted  the  current  versus  initial  angle  at  T = 1,  2, 

5,  and  10  shakes  for  altitudes  of  20,  30,  and  40  km,  with  correspond- 
ing curves  in  Figs.  11a,  lib,  and  11c  for  ionization  rate  (not  ioniza- 
tion). The  vertical  scales  are  arbitrary  but  compatible.  These  curves 
represent  the  contribution  to  current  and  ionization  rate  produced  by 
the  number  of  Compton  electrons  per  unit  angle  at  the  indicated  initial 
angle.  All  curves  peak  near  10  deg,  corresponding  to  the  peak  of  Fig.  6, 
but  the  varying  behavior  at  larger  angles  is  interesting.  At  the  highest 
altitude,  the  current  is  most  strongly  peaked  at  10  deg  at  early  times, 
whereas  the  peak  is  considerably  broader  as  time  increases  and  some  of 
the  electrons  which  started  at  large  angles  are  deflected  forward.  At 
lower  altitudes,  scattering  produces  a greater  broadening  of  the  dis- 
tribution. Energy  loss  effects  at  the  lowest  altitudes  cause  the  widest 
angle  (lowest  energy)  electrons  to  disappear,  as  can  be  seen  most  clearly 
in  Fig.  11a.  The  changing  shape  of  the  angular  dependence  as  time  in- 
creases makes  it  unlikely  that  it  would  be  a reasonable  approximation 
to  replace  the  Klein-Nishina  distribution  with  purely  forward-directed 
electrons,  as  has  been  done  by  other  codes. 

The  principal  results  of  the  computer  program  are  presented  in 
Figs.  12  and  13.  Figure  12  depicts  the  current  per  incident  gamma  ray 
(10  pcJ/2g  is  plotted);  Fig.  13  shows  the  ionization  in  thousands  of 
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secondary  electrons  per  incident  gamma  (10  Ng/g  is  plotted).  These 
are  shown  versus  time  at  altitudes  between  18  and  60  km.  Below  18  km, 
the  number  of  gammas  is  essentially  negligible  (see  Fig.  4).  Above  60 
km,  the  ionization  is  negligible  (Fig,  13),  and  the  current  is  virtually 
independent  of  altitude  (Fig.  12). 

At  high  altitudes,  the  effects  of  scattering  and  energy  loss  on 
the  current  are  slight.  Accordingly,  the  current  reduces  to  that  pro- 
duced by  magnetic  deflection  acting  on  the  Klein-Nishina  initial  dis- 
tribution. From  Fig.  12,  this  condition  essentially  holds  at  altitudes 
above  45  km  for  times  up  to  10  shakes.  As  the  altitude  decreases,  the 
•altering  begins  to  cause  the  electron  distribution  to  decohere.  The 
-a*  of  the  current  does  not  rise  as  high,  and  it  occurs  somewhat  earlier 
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Fig,  lib — Ionization  rate  v*  initial  angle; 
B=  0.6  Ey=  1.6  MeV  H = 30  km 
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Fig.  13  — Ionization  per  primary  gamma  ray;  B=0.6  Ey=  1 .6  MeV 
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in  time.  From  a peak  of  1.267  at  2.2  shakes  at  60  km,  it  drops  to 
1.255  at  2.0  shakes  at  40  km,  '1.06  at  1.6  shakes  at  30  km,  and  0.42  at 
0.8  shake  at  20  km.  The  peak  not  only  becomes  lower  and  earlier  but 
narrower,  with  the  decay  at  later  times  being  faster.  These  effects 
are  exactly  what  we  would  qualitatively  expect. 

A measure  of  the  effectiveness  of  the  scattering  in  reducing  the 
current  is  shown  in  Fig.  14.  Here,  we  have  plotted  the  ratio  of  the 
current  at  a given  altitude  to  the  current  at  high  altitude  (60  km). 

The  rapid  and  steady  decrease  of  the  relative  current  with  time  is  evi- 
dent. Thus,  at  20  km,  the  current  is  down  by  more  than  a factor  of  four 
for  times  later  than  2 shakes.  This  reduction  of  current  at  low  alti- 
tudes is  especially  important  because  the  fields  radiated  by  these 
currents  are  least  affected  by  attenuation.  The  initial  rise  of  the 
relative  current  at  high  altitudes  is  real.  It  corresponds  to  the  fact 
that  those  electrons  which  are  scattered  forward  spend  a greater  retarded 
time  radiating,  with  corresponding  field  enhancement,  than  those  which 
are  scattered  away  and  have  their  effective  radiation  reduced.  The 
effect  is  soon  washed  out  by  the  general  decoherence. 

The  ionization  versus  time  at  various  altitudes  is  displayed  in 
Fig.  13.  The  horizontal  line  at  23,600  represents  the  asymptotic  limit 
of  all  the  curves,  corresponding  to  complete  con\ ersion  of  all  the 
energy  of  the  Compton  electrons  to  secondary  ionization.  (The  mean 
energy  of  these  electrons  produced  by  the  first  scattering  of  the  1.6 
MeV  gamma  rays  from  the  bomb  is  800  KeV.)  Although  it  is  obvious  that 
the  electron  energy  is  lost  most  rapidly  at  low  altitudes,  because  the 
energy  loss  is  proportional  to  density,  the  curves  do  not  display  any 
striking  behavior,  except  to  note  the  large  range  over  which  the  ioniza- 
tion density  varies. 

If  we  refer  to  Eq.  (2.38),  we  see  that  the  ionization  density 
affects  the  fie1'1*5  by  it-a  appearance  in  conductivitv . where  it  is  divided 
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by  the  density.  Therefore,  we  have  plotted  10  N<,/pg  in  Fig.  15.  We 
call  this  quantity  the  equivalent  conductivity,  since  there  is  no  con- 
venient way  to  characterize  the  dependence  of  the  collision  frequency 
on  field  strength,  which  is  required  to  completely  specify  the  conduc- 
tivity. Figure  15  is  quite  remarkable.  At  low  altitudes,  the  change  in 
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the  Compton  electron  paths  due  to  scattering  and  energy  loss  causes 
the  secondary  electrons  to  be  produced  less  rapidly  In  retarded  time 
at  higher  densities.  Thus,  the  equivalent  conductivity  increases  with 
altitude  from  20  km  up  to  35  km.  Above  35  km,  the  ionization  lag 
(Eq.  (2.44))  becomes  the  dominant  effect,  and  the  equivalent  conduc- 
tivity decreases  as  the  altitude  increases  further.  The  combination 
of  these  effects  collapses  the  curves  of  Fig.  13  into  a very  limited 
region.  Thus,  if  the  curves  of  Fig.  15  were  replaced  by  a mean  curve, 
which  is  essentially  the  25  km  curve,  the  departure  from  this  mean  of 
all  the  curves  of  Fig.  15  for  all  altitudes  and  times  is  less  than  20 
percent.  It  accordingly  would  be  reasonable  to  replace  the  equivalent 
conductivity  per  electron,  regarded  as  a function  of  altitude,  by  an 
average  value  which  is,  of  course,  a function  of  time.  This  approxi- 
mation clearly  is  not  permissible  for  the  current. 

Figures  12  and  13  represent  the  current  and  ionization  per  gamma 
ray.  The  total  current  or  ionization  at  a point  is  found  by  multiplying 
these  curves  by  the  function  g(r)  (Eq.  (2.32)  in  general.  Fig.  4 for 
these  values  of  B and  E^) . We  shall  not  plot  these,  since  the  resultant 
curves  display  multiple  crossings  and  are  very  difficult  to  follow.  We 
note  that  the  curves  for  total  current  and  total  conductivity,  regarded 
as  a function  of  altitude,  peak  at  about  32  to  34  km,  corresponding  to 
the  peak  in  g(r). 

From  the  extensive  analysis  of  this  report,  we  draw  the  following 
conclusions : 

1.  The  development  of  the  electron  velocity  distribution  with  time 
has  been  determined  in  a consistent  and  satisfactory  manner 
that  includes  the  phenomena  of  magnetic  deflection,  energy 
loss,  atmospheric  scattering,  and  the  initial  Klein-Nishina 
velocity  distribution. 

2.  The  current  and  ionization  have  been  calculated  and  both  dis- 
play the  proper  qualitative  and  quantitative  physical  behavior. 
For  a delta  function  gamma-ray  source,  the  peak  of  the  current 
is  reduced  and  made  both  earlier  and  narrower  by  the  effects 

of  scattering,  whereas  the  ionization  is  determined  primarily 
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by  the  proportionality  of  energy  loss  to  density.  The  con- 
ductivity is  insensitive  to  altitude. 

3.  A computer  program  prepared  to  provide  numerical  values  for 
the  current  and  ionization  is  presented  in  the  appendix  to 
this  report.  The  program  operates  rapidly  compared  with  Monte 
Carlo  simulations,  thus  demonstrating  the  virtue  of  the  analytic 
approach.  If  the  program  is  augmented  by  an  integration,  arbi- 
trary gamma-source  time  dependence  can  be  included.  The  current 
and  ionization  as  calculated  here  can  be  used  as  inputs  for 
calculating  fields. 
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APPENDIX:  THE  COMPUTER  PROGRAM 


IMFLICIT  PEAL*9  (A-H,K,0-Z) 

DIMENSION  F (24)  , SEC  (111)  , P F:  (60) 

COMMON  FI  (SO,  121)  , ENN  (1  1 1)  , AJJ  (1  1 1)  , PP  (6,121  ) 

R,F (200)  ,CC(121)  ,SS(121),TT  (5,  121)  ,DP<5,  121)  ,00(50) 

6,  ENN.MP  (3,  111)  , A J J M K (3,111)  , l\'  1 (5,  121)  , EN2  (5,121)  , EK  3 (5,  121) 

6, ENNM  (3)  , AJJM  (3)  , AJ  1 (5, 121 ) , A.T2  (5,  1 21 ) , AJ 3 (5,  1 2 1 ) 

6,  ENK  (3)  ,AJK  (3)  ,SNA  (3,9) , AJA  (3,9) 

6,DEA(9)  , Q I (9)  ,v»C(Q)  ,Jfi(9)  ,QS  (9)  , PU  ( 1 2 1)  , DEQ  (5)  , DEB  (5) 

6,CK  (121)  , FP  (14,  121)  ,GI(5)  ,WI(5>  , DE(5,4)  ,AK1  (5,4) 

£,AK2(5,4)  , AP 3 (5,  4)  , E"  2 (5,4)  , F .M  3 (5 ,4 ) , Ev  M ( 3)  , EX2  (4)  , EX  3(4) 

E,FNP  (1 1 1)  , AJP  (1  1 1) 

C-iHICN  I!i),DELH,YY,YYnEL,VV,H,2,GG,KK,ni,FMY,X,XSN,FHH,B,C,YO 

£,TTP,TTO,TTI,TTU,TTV 

6,  MLIM,  MLIMK,  Ih,I.V,I  Y«L,IYH  , IYL.LGO,  IYV,  MVX 
G(ZDIIM)  = GG*i)EXP  (-KK*Z  DUP) 

UO  (ZrUB)  = ((G(ZDUM)  -IDO)  *(GG*1D'')/  ( (G  <ZTMIM)  ♦100)  * (GG-1D0)  ) ) 

1 **2.164D0 

A (ZPU M)  =2. 164D0*FK*G (ZDUV ) / («  ( Z DM M ) ** 2- 1 DO ) 

0 ( Z DUB)  = 5.6T>569Dt*DS.lRT(G  (ZD'1.1)  **2-U0)  *UU  ( ZDOfl ) /BB 
ELL.  (MDM)  =.  3 45  34  903720  9 3D'  * (IDO*  . 16D0*DI  OG  (IDO  *3  . 3 1D4*HHH*B*NDN)  ) 
DATA  P/.  ICC  3DO,  .072  54  DO,.  052  6 5DO,  .0  3-33  ICC,  . J2  79  7P0.  . C 204  6»r  0 
6,  . 0150  3D  0,.  01 106D0,  . e»  1 S 3 L -2, . 6 024  P-2,  .44661  -2,  . 3323D-2,  .2486D-2 
6,  . 187 ID- 2, . 14190-2, . 1 0840-2, . 8 34  8D-  3, .6481  P-3 , . SC66D-3, . 399  3D- 3 
6,  . 3142D-3,  .24-33D-3,  . 1 962!  -3,  . 15441-  3/ 

DATA  ZEHO,P,u/0DG, 3. 1 4 1 69 2 65 3 59979D0 , 1 . 74 5 32925 1 994 33 D- 2/ 

0, 0*1, ON/0.9510  56  51629515D. 1,0.5677 -3  5252292470  0/ 

FPn  3=4  DO/ 3DC 
DCON-1 . 7435L4/24Tr' 

c *****  BF  = MAGNETIC  FIELD,  E E - GAP*' A ENERGY, 

C *****  5 = TIMF-STEP,  C = ANGLE-STEP 

FSAD  1, BB,FE,B,C,Y1T 
1 Y MI  = YM I 
1 FOP  NAT ( 6D  12 . 4) 

MLI  N=  1 1 DO/  B*  1 . C 1 TO 
*L  I MK  = P!  LI v - 2 
K=EE/0. 51100400 
K 1 = K ♦ 1 D 0 
A2=K*2 DO 


KS=F*K 

K1S=K1*K1 

K2S=K2*K2 

S= (2D0+K*  (HCO*K*  (9D0  + K) ) ) / ( K * ( 1 DO *2 DC*K ) ) **2 
1 - (1D0*K*  (1D0-r. 5 DC  *K) ) *DLOG ( 1D0+2D0*K) /F**3 

PELH=2DC 

READ  1,HH0,YY0,III!1,  YY1 
IIH  = HHO 
I li  = H I1 1 

IP  (IH.NE.O)  II1  = 1H-1 
10  I H = I H ♦ 1 


r H H =r  (IH) 

KP=6.5954D0*PHH/DP 
r^INT  11,HH.BB,EF,6,C 
11  FOF  n at  ( 1 Si  1 , 1 ALTITUDE  =• 

6 GAMMA  ENERGY  = • , F6 . 2 , /2 1 X , 
C * , F6 . 2, * DLG . 1 ) 

LO  512  1=1, MLIM 
AJJ  (PI)  = ZFEO 
ENN  ("*)  = ZERO 
ENP (M) = Z EPO 
AJP (M) =ZEHO 


,F4 . 0,  • K.N 
• TIKE-STFE 


.,  MAGNETIC  FI  El  D =',F6.2,', 

= ’ , P(; . 2 , ' SEC.,  ANGLE-STEP  = 


^ RECEDING  PAGE  BIaNK-NOT  FILMED 
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512  CONTINUE 
TTV  = 1 5D0 
HVX=151 
Y YDEL=3D0 
YY= YYDEL 

20  CONTINUE 

c*******  PART  2 

EHO=ZERO 
ZZ=YY*Q 
XSN=DSIN(ZZ) 

X=DCOS (ZZ) 

XS=X*X 

XSN  = 1 DO-X  S 

XS3H=1D0-3D0*XS 

XS4N=1D0-5D0*XS 

XS5H=3D0-5D0*XS 

01(1)  =XSN*Qf' 

QI  (2) =XSN*QN 
QI (3) =Z  ERO 
01  (4)=-QI  (2) 

0I(5)=-QI(1) 

KX  = (K*X/K 1) **2 
XSK  2=K*  XS  *K  2 

(1D0»KX)/(1D0-KX) 

SS(1) =Z  ERO 
CC  ( 1)  =ZEPO 

BU(1)=DS0RT (IDO- IDO/ (33*33) ) 

OK ( 1) = Z ERO 
DO  21  1=1,5 
TT  (I, 1)  =Z  ERO 
DD (I, 1) =1 D0-X*BU  (1) 

A J 1 (I, 1) =ZERO 
A J 2 (1,1) = ZE*0 
A J 3 (1,1) = ZF~0 
EN2  (1,1)  = ZERO 
EN 3 (1,1) =ZEfO 
0II=QI(I) 

31  (I) =01 1/XSN 

41(1) =2D0/(1D0*GI(I)  ) 

0I2=QII*0II 
Q8 (I) =1 D0-5D0*QI2 
QS  (I)  =1  DO- 3 D0*Q  12 

21  OC (I) =QII*  (2D0+QR (I) ) 

EHN2=(K1S-XSK2)/(K1S-KS*XS) 

EKH1  = EN!12*1C0/E!1M2-  (1D0-XS)  * ( 2 DO* K 1 * X/ ( K 1 S- XSK 2)  ) ** 2 
ENY=2D0*K1S*X*DSyRT (1D0-XS) *Efitf1/S/(K1S-XS*KS)  **2 
rRINT  1 2 , Y Y , EPY 

12  FORHAT(//'  INITIAL  ANGLE  = • , F4 , 0 , ' DEG.,  R5L  INITIAL  NUHBER  =• 
1,F9.5) 

V V = DLOG (GG) /KK/0 
VV=VV-C 

IF  (VV.3T. 119D0)  VV= 1 1 9D0 

ITS  X=TV/C ♦ 1 DO 

IYHXH=I YNX-1 

IF (IYHX . LE. 5)  GO  TO  96 

YO=  ZEPO 

H=C 

«=?V 

IYNL=IYNX 
L»  = 2 


r 
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LGo  = 1 

C ASSIGN  6 02  10  LGO 

CAI.I,  PA  ST 7 
502  CONTINUF 
Y = ". 

I Y = 2 

25  CONTINUE 

CC  (IY)  =FI  (1  , IY) 

SS(IY)  = EI(2,IY) 

Z = Y*0 
GOZ=G  (2) 

OK  (IY)  = . 36207D0*  (G3-30Z  ) /FHli  -X*CC(IY) 

BU  ( IY)  = DSUrT  ( 1BC-1  DO/C.nz  **2)*UU(Z) 

Y-Y  *H 
IY=IY*1 

IF  ( Y . LE . V V)  GO  TO  25 

I YU  =IY 

IYV=IYU-1 

IY2=TY- 2 

I Y 3 = 1 Y-  3 

iyu=iy-u 

I Y 5 = IY- 5 

OK  (I YU)  =0K  (IY5)  -5nO*(OK  (j.Y4)  -OK  (IYV)  ) ♦1000'-  (OK  (IY3)  -OK  (IY2)  ) 

BU  (IYU)  = BU  (IY5)  - 5D0*  (BU  (IY4)  - EU  (I  YV)  ) +10nC*  (UU  (IY3)  -n'l  (IY2)  ) 

SS  (IYU)  = SS  (IY5)  -5D0*  (SS  (IY4)  -S3  (IYV)  ) ♦ IOd''*  (SS  (IY3)  -SS(IY2)  ) 

CC  (I  YU)  = CC  (IY5) -5D0* (CC  (1 Y4) -CC  (I YV) ) *1000* (CC (IY3) -CC (I Y2)  ) 

TTP  =0K  (IYU) 

TTO  = XSN*SS  (IYU) 

TTI.  = TTP-TTO 
TTU  =TTP  *TTO 


L V = 50 
LGO  = 3 

C ASSIGN  640  TO  LGO 

CALL  PAPT7 
540  CONTINUE 

C*******  EA?T  3 

DO  3 C 1 J = 1,  121 
DO  301  1=1,5 
30  1 PP  (I,  J)  = Zt’B0 
DO  302  1=1,14 
302  FF  (1,1)  =ZZFO 
Y = H 
IY  = 2 
I YL  = 1 20 

IH=QMIY)  ♦ X SN*S3  (IY) 
34  CONTINUE 
TG  = TH 
Z = Y*0 

SI N Z = DS I N (Z) 
COSZ=t)CCS  (Z) 

1=1 

36  CONTINUE 
CCI =CC ( I Y ) 

SSI=SS(IY) 

E 1 1 Y = EI  ( I , I Y) 
EI1IY=EI (1*1, IY) 

LI2 I Y=EI (1*2, IY) 

El  3 IY  = EI (1*3, IY) 

EI4 I Y = EI  (1*4, IY) 


oe>  c*  o c*  c*  ©»  cj  o C-  c*> 
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F(I)=EIIY*CCI-EI1IY 
FI  = F (I) 

F (I ♦ 1 ) =EIIY*SSI-EI2IY 
FI  1 =F  |IH) 

F (I+2)=SSI*FI-CCI*EI2IY*EI3TY 
F (I *3) =CCI*  (FI- Ell I Y)  *EI4 I Y 
F ( I ♦ 4 ) = SSI* (FI1-EI2IY) + EI (1+5, IY) 

F (I *5) =CCI* (CCI* (FI-2  D0*E 1 1IY ) »3D0*EI4IY) -El (I+6.IY) 

F ( I *6) =CCI*F(I*4) -SSI*  (SSI*EI II Y- 2D0 *EI3I Y)  -EI(I*7,IY) 

F (1*7)  = SS I * F (I ♦ 3 ) -CCI * (CCI*EI2IY-2D0*EI3IY) -El (I*8,IY) 

P ( I *8)  = SS  I*  <iJS I*  (FI  1-  2D0* EI2I Y)  ODO* El  ( 1+  5,  I Y ))  -El  (I*  9,  I Y) 

F (I *9) = ZERO 
I=I*10 

I F ( I . LT . 42)  GO  TO  36 
PP 1 = F ( 1 ) 

PP ( 1, IY ) =PP1 
PP2-F  (2) 

PP ( 2, IY ) =PP2 
PP  3 = F ( 3) 

PP ( 3, IY) =PP  3 
P P 4 = F (4) 

PP (4,IY) =PP4 
PP5  = F (5 ) 

PP(5,IY) = PPb 
DO  31  1=33,49 
31  F(I*20) =F  (I) -F  (1-20) 

Fv (1 ,IY)  = (2 DO*  (F  (14) ♦ F (1b) ) +XS3M*PP4) /3D0 
FF ( 2, I Y ) =2DC*X*PF  3 
FF ( 3, IY ) =PPb/3D0 

FF(4,IY) =X*  (F  (31) * (20  0*PP4-PI 5)  *3D0*F(32) *PF3-F  (56) -F  (57) 

6 ♦XS5«*(F(41)  *PP4-.5D0*F(66)  )) 

FF (5,IY) =-3D0*F(31)  *FP3-F (3  2) * ( 2D0*  PP5- PP4) +F  (5  8) + F (59) 

& +XS4M* (-F (41)*PF3-F(4  2)  *PP4*1 .5d0*F (h8)  ) 

FF  (6,IY)  = X*  (F  (4  1)  *PF5*2D0*F  (4 2)  * PP3-  1 . 5 DO*  F (6 7)  ) 

FF (7,IY) =-F (42) * PF5  + 0 • 5D0*F  (69) 

FF(8,IY)  =2D0*(-SINZ* (2D0*F (11) ♦ XS  3F *PP1 ) ♦2D9*COSZ*F  (12) )/3D0 
FP (9,IY)  =2DC*X* (-SINZ*FP2+COSZ*PF1) 

FF ( 1 0,  I Y) =2D0*COSZ*PP2/3D0 

FF  ( 11,1  Y)  =0.2D0*X*(SINZ*  (2 DO*  ( F ( 1 4)  ♦ 2D0  * F (1  5)  ) ♦XS5!1*F  (24)  ) 
♦2DO*COSZ*F  (13) ) 

-(0.8D0*SINZ*X* ( 4 D 0 * F ( 3 1)  *PP1 ♦ 3 DO *F  ( 32)  *PP2-2D0*F  (54) 

-1 . 5D0*F  (55) ♦XS5«* (2D0*F  (41) *PP  1-F  (64) ) ) 

♦ 0.  8D0*COSZ * X* (2D0*F( 31) *PP2-3D0*F (32) *Pr 1>  0. 5D0*P  (53) ) ) 

FF (12, IY) =0.4 DO* (SI HZ*  (F(13) -XS4H*F  (23) ) *CCSZ* ( F ( 15) ♦2D0*F(14) 

♦0. 5D0*XS4H*F (24) ) ) 

- (0.8D0*SINZ*(-3D0*F( 31) *PP2* 2DO *F ( 32 ) *PP 1 ♦ 0. 5D0*F (53) 
-XS4M*  2D0*(F(41) *PP2*P  (42) *PP1  -F (63)  ) ) 

♦ 0. 8DO*COSZ* (3D0*P(31) *PF  1 + 4D0*F  (32) *rF2-2D0*F (55) 

-1. 5D0*F(54) ♦XS4F* (2D0*F(41)*PP1-F  (64) ) ) ) 

FF ( 1 3,1  Y)  =0 . 2D0 * X* (SINZ*F (25) - 2 DO*COSZ* F ( 2 3) ) 

- (0.8D0*X* (SIKZ* (2D0*F  (42) *PP2-F  (65)) 

-2D0 *COSZ * (F (41) *PP2+F (42) *PP 1- F (63 ) ) ) ) 

PF  (14, IY) =0 . 2DO*COSZ* F (25) 

6 - (0.8D0*COSZ*(2D0*F(42) *PP2-F (6  5)  ) ) 

Y = Y*H 
I Y= I Y* 1 

TH  = QK(IY)  +XSN*SS  (IY) 

IF(TG.LT.TTL.XND.TH.GT.TTL)  1YL=IY 
IF(t.LE.VV)  GO  TO  34 
PRIBT  853,IYL*IY0,TTI.  ,TTU  , TTP, TTQ 
853  FORMT</'  I YL  =•  , 14  , • , IYU  s',14,'. 


TTl  =• , F9.«, ' , TTO  *• 
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6F9.4,',  TTP  =',F9.4,»,  TTO  =',F9.4) 

L V = 42 
LGO  = 2 

C ASSIGN  632  TO  LGO 

CALL  PART7 
532  CONTINUE 

C*******  PART  4 

DO  401  1=1,200 
401  F ( I ) =ZF 10 
I = H 

iy=2 

70  CONTINUE 
Z=Y*Q 

SIN Z = DS IN  (Z) 

COSZ  = DCOS  (Z) 

1=1 

43  CONTINUE 
SSI=SS(IV) 

CCI =CC (IV) 

El  I Y = EI  (I  , I Y) 

Ell IY  = EI  (1*1, IY) 

EI2  IY  = EI  ( 1 + 2, IY) 

EI3IY=CCI*EIIY 
EI4 IY=SSI*EIIY 
EI7 IY=EI  (1*6, IY) 

EI5IY=CCI*EI71Y 
EI6IY=SSI*EI7IY 
EI8 IY  =F I (1*7, IY) 

EI9IY=EI (1*8, IY) 

F ( I ) =EI3IY-F.11IY 
F (I ♦ 1 ) =EI4IY-FI2IY 

F ( I ♦ 2) = CC I*  F (1*1) -SSI*FI1 IY+FI  (1  + 3, IY) 

F (1*3)  =CCI*  ( F ( I)  -El  1 1 Y ) *FI  (1*4,  IY) 

F (I*4)=SS1*  (F  (I*  1)  -EI2IY)  *EI  (I*5,IY) 

F (1*5)  =FI5IY-EI81Y 
F <!♦*»)  =EI6IY-EI9IY 

F (I *7) =CCI*F  (I ♦6) -SSI*EI8IY*EI  (I*  9,IY) 

F (148) = CCI*  (F  (1*5) -EI8IY)  *EI (1*10, IY) 

F (I *9) = SSI*  (F  (1*6) -EI9IY)  ♦ El  (I*  11 ,IY) 

F.IC1Y  = E1  (I*  1 2 , I Y) 

F (I *10)  =CCI ♦"ICIY-EI  (1*13 , I Y) 

F(I*11)=SSI*EICIY-EI(I»14,IY) 

E1FIY  = EI  (1*15, IY) 

F (I  ♦ 1 2)  =CCI * EIFI Y- E I (I*16,IY) 

F (1  + 13)  =SSI*EIFIY-EI(I*17,IY) 

E I 1 HI Y=  El  (1*18, I Y ) 

F (I *14)  =CCT*FT1RIY-F1  (I*19,IY) 

F (1*15) =SSI*ET1RIY-FI  ( T*2  C , I Y) 

1=1*21 

I F ( I . LT  . 2 3)  GO  Tf  4 3 

FF  (4,  IV  ) = (-  (X*  (4l'0*F(4)  *2P0*F  ( 13)  *3D0*  (F(5)  *F  (12)  ) ♦ F ( 8)  - F < 1 5) 

6 *XS  5F  * ( 2D0*  F (25)*F(34))))  *FF(4,IY))*O.8D0 
PF  (5,IY)  = (3D*  (F  (9)  *F  (11)  ) *F  (3)  *4D0*F  (10)  - F (14)  *2J)0*F  (1b) 
f,  *XS4*,*(4D0*F  (24)  *2D0*  (F  (3  0)  *F  (32) ) *F  (35))  *FF  (5,IY) ) *C.8D0 
FF  (6,  IY)  = (-  (X*  (2TO*  (F  (26)  *F  (3 3)  ) *4D0*P(2°)  ♦ F ( 36 ) ) ) ♦ FF  (6 , 1 Y)  ) * . 8DP 
FP  (7,IY)  = (2 10* F (31)  *F  (37)  *FF(7,IY)  ) *C.fiD0 
FP(  11,  IT)  =-  (X*  (-SINZ*  (4D0*F  (1)  *3DC*F  (7)  *2  DC *XSS P*F  (22) ) 

6 ♦COSZ*(3D?*F(2> -2DO*F  (6 ) ) ) ) *0 . 8 D0*F F ( 1 1 , 1 Y) 

FP(12,IY) =-  (SINZ*(3D9*F(6)  -2 PC*  (F  (2) *XS 4fl * ( F (23 ) *F (27) ) ) ) 
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6 -CO  SZ  * ( 4D0*  P (7) *3D0*F (1) ♦2D0*XS4K*F  (22) ) ) *0 . 8D0 *FF ( 1 2 ,1 Y) 

FF (13.IY) =-  (2D0*X* (-SINZ*F(28) YCOSZ*  (F (23) *F (27) ) ) ) *.8D0*FF  (13,IY) 
FF(14,IY)=  2DO*COSZ*F  (28)  *C . 8DC+ FF < 14, I Y) 

Y = Y 

IY=IY*1 

XF(Y.LE.VV)  GO  TC  70 
DO  411  1=1,  14 

411  FF  ( I, IY ) =FP(I,IY5)  -5  DO* ( FF  (I,IY4)-FF(I,XYV)) 

6 ♦10D0*(FF(I,IY3)-Ff(I,IY2)) 

c *******  rAFT  4A 

I Y = 2 

Y = H 

26  Bn  I Y = BU  ( I Y) 

Z = Y *0 

SINZ=DSI»  (Z) 

COSZ  = DCOS  (Z) 

QKIY=0K  (IY) 

&L  = SS  (I  Y) 

DO  49  1 = 1,5 
011=01(1) 

0RI=QR(I) 

0SI=0S(I) 

gci=QC(i) 

TT (I,IY) =QKIY  ♦0L*QII 

A.T1  (I,  IY)  = <X*SI«TZ*QII*COSZ)  *BUIY 

DC  ( I , XY) =lrO-BUIY  * (X*COSZ-QII*SINZ) 

EN2  (I,IY)  =FF  (1 ,IY) ♦ QI I*FP  ( 2 , I Y)  ♦QSI*FF  (3, IY) 

EN3 (I,IY) =FF (4,IY) ♦QII*PF (5,IY) +yRI*FF(f , IY) ♦QCI*FF (7,IY) 

AJ2 (I,IY) =(FP(8, IY)  ♦0II*FF  (9,IY) +0SI*FF (10,IY) ) *BUIY 
4 9 AJ3  (I,IY)  = (FF  (1 1 , T Y ) *QII*FF  ( 1 2 , I Y) ♦ QHI *FF ( 1 3, I Y ) +QC I*PF  ( 1 4, I Y) ) 

6 *BU I Y 

Y = Y*H 
IY=IY»1 

IF  (IY.LE. IYU)  GO  TO  26 
CALL  PART5 

IF  (YY.GE.29.9D0)  Y Y DE  L=6D 0 
YY= YY+Y YDEL 

IF(YY.LT.89.9D0)  GO  TO  2C 
PRINT  912, YY 

912  FORMAT (//*  INITIAL  ANGLE  =',F4.C,'  DEG.1) 

96  CONTINUE 
".=  1 

T=ZEFO 

906  CONTINUE 

TVT1  = . 5D0*  (TTV-T)  / TTV 

ENN  (B)  = ENN(S)  ♦ TVIJ*ENP(P)*FP03 

AJ J (H)  = A JJ  ( H)  ♦TVU*AJP  (B) 

IF  ( (2* ( (d-1)/2) . EQ.B-1 .AND.T. LE.2.01D0)  .op. 

6 (5*(<B-1)/5) .E0.B-1.AND.T.GT.2.01D0)) 

CPKINT  855,T,AJJ  ( fi)  , ENN  (H) 

855  F0PHAT(F8.2,2(37X,  F10.5)) 

T = B*B 
B=B  ♦ 1 

IF(B.LE.BVX)  GO  TO  906 
574  FOFBAT (1H  ) 


C******* 


PART  6 
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66  1 
03 


05 


1 C 1 
102 
1 C 3 

104 

105 

106 

110 

120 

664 


D = DC0N*«‘HH*  e 
P^INT  661 

FOFI*AT(/PX,  'TIME*  ,8X,  TRIDEL*  ,0X,  'SECOEL'  ,9X,  • JPUIDEL') 
« = 1 

CONTINUE 
T = M*B 


M = " ♦ 1 

M A X = (1-  1)  /2  ♦ 1 
DO  95  1=1, TAX 
ELLI=FLL  (1-1) 

I P ( ELLI  . GT . IDO) 
ELI. N = £LL  (1-1) 

IF (ELLM .GT.  IPO) 
c ? ( I)  =E  LL  M*r  NN  ( 
CO*'TlNU  E 

IF  (1. or. 7)  00  TO 
I F ( X . O.T  . 6 ) GO  TO 
IF  (H.GT.5)  00  TO 

I F ( 1 . GT  . 4 ) 00  T"' 

IF (M.GT. 3)  00  TO 

IFJ1.GT.  2)  C-C  Tf1 
cNNJ1=12D0*Fr(1) 
GO  T°  120 


ELL  I = 1 DO 

ELL  1= 1 DO 
I)  +FLLI*  FNN  ( 

106 
105 
104 
103 
102 
10  1 


to- I ♦ 1 ) 


?NNJ1  = 8D0*',r  ( 1 ) ♦ 16P0*EP  (2) 

GO  TO  12P 

I NNJ1=9PC*F’-  (1)  ♦27D0*PF  (2) 

GO  TO  120 

FNN«1H=8PC*  ('  * (1)  *410 *7F  (2)  *F  ■ ( i)  ) 

GO  TO  120 

ENNJM=8D0*T I (1) ♦ 31  DO* PE  (2)  +21  CO*?”  (3) 

GO  TO  120 

t NN J1=8  DO  *l  : ( 1 ) ♦ 3 1 DC*  FF  (2 ) ♦ 2 0 D 0*0  c ( 3)  ♦ 1 3DP*  I-  c (4  ) 
GO  TO  1 2P 

ENNJ1=8D0*’,r  ( 1)  ♦ jl  r>P*PP  (?)  +2PD0*0?  ( 3)  ♦2  5DC*F'  (4) 
IE  (1.Ew.°)  GO  TO  120 
LO  11n  j = 5 , 1 A X 
ENNJM=ENNJ1*24P0*?F  (I) 


C o ’)  T I *’  U F 

IF  (2*  (1/2)  . HE.  M)  ENNJH="N  NJW-1  2DO*FF  (1,  AX) 
CONTI  V II  E 


SEC  (1)  = D*  EN  "M  M 

OF  I NT  6 6 4 , T , E N N ( M)  ,SEC(M)  , A.KI  ( I- ) 

FOR1AT(F12.2,F15.5,F15.2,2F15.5) 

IF  (fl.  L-".  II  1 1 K ) GO  TO  93 


1!H  = HH*DF.LH 

IF (HH.LE.bOPP)  Gf  TO  10 

CALL  EXIT 

END 


SUBROUTINE  taFT5 

c, ******  PART  5 PEVTSEP 

IMPLICIT  FEAL*0  (A-H,K,0-Z) 

CO  IKON  El  (50,  121 ) , ENN  (111)  , A.IJ  ( 1 1 1 ) , PP  ( 5 , 1 2 1 ) 

&,F  (200)  ,CC(121) ,SS(121),TT(5,  121)  ,PP(5,121)  ,Q0(50) 
6,ENN(1K(3,111)  ,AJJf1K(3,111)  ,EK1(5,121)  ,EN2(5,121)  ,EN3(5,121)' 
6,F.NVM(3)  , AJJE  (3)  ,AJ1(5,121)  ,AJ2(5,121)  ,AJ3  (5,121) 
fi,  ENK  (3)  , AJK  ( 3)  ,FNA  (3,9)  , AJA  (3,9) 


r 
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5. DDA  (9)  , QI  (9)  ,QC  (9)  ,CH(9)  ,0S(9)  , BU  { 1 2 1)  , D EQ  (5)  , DEF  (5) 

E.QK  (121)  , F»  (14,  121)  ,GI(5)  , WI  ( 5)  , DE (5, 4)  ,AK1  (5,4) 

6,  AK2  (5,4)  ,AK3(5,4)  , PH2  (5,4)  , EM 3 (5,4)  ,ENH(3)  ,EX2  (4)  , PX3(4) 

6 , EN  P ( 1 1 1)  , AJP(III) 

CON  HON  HH,DELH,YY,  YYDEL,VV,1I,N,GG,KK,  BP , EH Y , X , X SN ,F HH , B , C, YO 
&,TTP,TTQ,TT1,TTU,TTV 

6,HLIH,NLIHK,IH,LV,IYHL,IYU, IYL, L30, IYV, HVX 
G(ZDUH) =GG*DEXP(-KK*ZDUH) 

UU(ZDUH)=  ((3(ZDUH>  - IDO)  * (33  *1  DO)  / ( ( 3 (ZDUH)  ♦IDO)  * (33- IDO)  ) ) 

1 **2 . 1 64  DO 

A (ZDUH)  =2.164D0*KK*G(ZDUE)/(G  (ZDU") ** 2- 1 DO ) 

O(ZDUH)  =5. 68569D0*DS0PT  (3  (ZDUH)  **2-1  DC)  *UU  (ZDUH)  /BB 
DATA  ZERO,P,Q/0D0, 3. 1 4 1 59 2653 5 897 9D0 , 1 . 745 32925 1 994 1 3 D- 2/ 

&, QN, ON/ 0.9510 56 51629515DC, 0.5 87 785252 29 24 7d0/ 

FRO 3=4D0/3u0 
PFINT  551 

551  FOp  H AT (/4  X, 'TIHE',4X, • JPH 1 1 • , 4X , • JPHI2 • , 4 X, • J PH  II • , 5X , ' JPHI • , 4 X , 
6'CUHJPHI* , 5X, ■ PF11 »,5X, *PRI2 ,,5X,,PRI3*,6X, »PRI' ,5X, *CU 

&HP7I ' ) 

Dn  511  1=1,200 
511  F ( I ) =ZE SO 

DO  513  1=1,3 
ENNHK (1,1) = ZERO 
AJJHK  (1,1)  =ZFPO 
ENNH  (I)  =ZES" 

AJJH  (I)  =Z  F70 
513  CONTINUE 

facqh=ehy*yydfl*q 

ttub=ttu-b 

TTUB2  = TTU  B- P 
TTLB=?TL-B 
TTLB2=TTLB- P 
BDN  1= 12D0  *B 
EDN  2 = BDN 1 *B 
BDN  3 = 2D0*  B* * 3 
HLU=TTL/B> 1 .01  DO 
NLX=NtIH 

IP(HLU. LT.HLIH)  HLX  = HLU 
HU  U=TTU/B ♦ 1 . C 1 DO 
HUX=HLI HK 

IF(HUU. LT.HLIHK)  HUX=HUU 
H = 1 

80  CONTINUE 
T = N*B 
H = * ♦ 1 

DO  81  J = 1 ,3 
ENF (J) =ZERO 
8 1 A JK  (J) = ZE  SO 
1=1 

82  CONTINUE 
IT  = 0 

84  IY=IY*1 

IP (T.GT . TT (I, IY) ) GO  TO  84 
SLD= . 2D0 

I P ( I Y . GT . 2)  I T=I  Y-  1 
IP  (IY.3E.IYV)  IY  =IY -1 
TTA  =TT (I,IY-1) 

TTB=TT (I , IY) 

TTC=TT (I, IY*1) 

TTD=TT(I,IY*2) 

TTE= (-TTA*9D0* (TTB*TTC) -TTD) /16D0 
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TTF  = (TT A-27P0*  ( TTB-TTC)  -TTP)  /24D0 
TTG  = (TTA-TTF-TTC+TTD)  /4DP 
TTH  = (T-TTE)  /TTF 

TTI  = (-TTA+3  CO*  (^TP-TTC)  + TTP)  /tiuO 
TTJ  =TTG/TTF 

rTK=2P0*TTJ*TTJ-TTI/TTF 

AP  = . SO0  ♦TTF*  (IDO-TTIi*  (TTJ - T TK + TTII)  ) 

A PM  2 = AF "2  CO 

AP5M=AP*AP-1D0 

AL1=-AP* (AT -IDO) *A^M2/6nO 

AL2=0.5nG*ArSM*APM2 

AL  l=-0.  5D0*AP*  (A  P+1  DC)  *AT  M2 

AL4=AP* APSM/6DO 

DDA(I)=AL1*DP(I,IY-1)+AL2*DI'(I,IY)+AL3*PC  (l,IY+1)+AL4*PD(I,IY+2) 
ENA  (1,1)  = IDO 

SNA  (2,1)  = AL  1*5N2  (I,  I Y - 1)  ♦ AL2*E"2  (1,IY)  ♦ Al ’*?::; 2 ( I»IY*1 ) 

£ +AL4*LN2  (I,  TY  + 2) 

ENA  (1,1)  =AL1*EN3  (I,  IY-1)  +AL2*EN3(I,IY)  ♦ AL  3 *?N  1 ( I , IT  ♦ 1 ) 

£ +AL4*EN3  (I, IY  + 2) 

AJA  (1,1)  =AL  1*AJ1  (I,  IY-1)  ♦ AL2*AJ  1 (I  , I Y)  ♦ AL3*AJ  1 (I.IY+1) 

6 +AL4+AJ1  (I, IY+2) 

AJA (2,1) = AL  1*AJ2  (I, IY  -1)  ♦ AL2*  AJ 2 (I,  IV) +AL3+AJ2  ( I , IY+ 1 ) 

£ ♦ AL  4 * AJ  2 ( I,  I Y + 2 ) 

AJA (3,1) = AL1*AJ3 (I, IY-1) ♦ AL2*AJ3  ( I , I Y ) +AL3*AJ  3 (I, IY  + 1) 

£ +AL4* AJ  1 (I, TY+2) 

.1  = 1 

PP A I = DD  A ( I) 

66  ENK  (J)  = FNK  (J)  ♦ SID  *ENA(J,I)/DDAT 
A J K ( J)  = A J K ( J)  ♦ SJH  *AJA  (J,  I)/DPAI 
J =J  ♦ 1 

IF (J.LS.  1)  JG 
1 = 1 + 1 

IF(T.LE.S)  ( 

DO  47  J=1 , 3 

AJJMK  (J,“)  ' I)  C*  F ACQF  * A JK  (J) 

P 7 F»|NMK(J,M)  = VAC*.3*r.NK  (J) 

IF(N.LT.MLX)  GO  ”0  80 

90  C O N T I *1 0 E 
T=  Z Fr0 

M = 1 

ENNM  (1)  = F AC  C'N/  ( 1 rO-X*BU  (1 ) ) 
r,r  Tr'  8 5 
9 1 M = 2 

r=? 

CNN"  (1)  =ENNMK  (1,2) 

AJJ*  (1)  = A J J * * (1,2) 

ENNM(2)  = (in0*ENFMK(2,1)-20P0*FNNMK(2,2)+4D3*FNNMK(2,3) 

6 + 4DC  * EN  'i  M K (2,4)  -ENN""'  (2,5)  ) /BDN2 

F NN  M ( 3)  = - ( 3 DO*  ENNMK  (1,1)-10C0*ENNMK  ( 3, 2)  ♦ 1 2D0  *ENNMK  (3 , 3) 

F -6D0*ENNMK  (3,4)  + ENNMN  (3,5)  ) /BDN  3 
A.TJ"  (2)  =-  (3P0*AJJMK  (2,  1)  ♦ 10PC*AJJMK  (2,2)  -18D0*AJJNK  (2,3) 

£ +6  P 0 * A J J M K (2,4)  -AJJMN  (2,5) )/BEN1 

AJJM  (3)  = (1 1 DO*AJJMK  (3, 1) -20DC*  AJJMK  (3,2) +6D0*AJ JMK  (3, 3) 

£ +4D0*AJJMK  (3,4) -AJJMF ( 3,5) ) /BDN 2 
GO  TO  85 
92  T = M*B 
M = M + 1 

ENNr.(l)  = ENN  *K  (1,*) 

A JJ  M (1)  = A JJ  MK  ( 1 , M) 

I F (T . GT  . TTL62)  Of’  TC  901 
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ENNH(2)  = (-ENNHK (2, N-2)  ♦ 16  DO  * { ENNHK  (2 , H- 1)  ♦ENNMK  (2,fl*1)) 

6 -30D0*ENNHK  (2,  ) - ENNMK  (2,  M+2)  ) /BDN2 
ENNH  (3)  = - (ENN"K  (3,  H-2)  -2D0*  (ENNHK  (3,M-1)  -ENNHK  (3,  (1*1)  ) 

6 -ENNHK(3,H*2))/BDN3 

AJJH  (2)  = - (-AJJHK (2. H-2) *8D0* (AJJHK (2, H-1)  - AJJHK  (2,M*1) ) 

& ♦ AJ  J HK  (2,H*2))/BDN1 

AJJH  (3)  = (-AJJHK  ( 3,  H-2)  ♦ 1 6 DO  * ( A J JH  K ( 3 , M-  1)  ♦ A JJ  (IK  ( 3 , H ♦ 1 ) ) 

6 -30D0*AJJHK (3,H)  -AJJHK  (3,  ft*  2)  ) /BDN2 

50  TO  85 

901  I F (T . GT . TTL  P ) GO  TC  902 

FNNH(2)  = (-F.NNHK(2,H-3)  *4D0 ‘ENNHK  ( 2,  H-2)  *6D0*ENNHK  (2,H-1) 

E -20D0*EHNMK  (2,H)  ♦ 1 IDO* ENNHK  (2,  H*1) ) /BDN2 
FN  N fl  ( 1)  = (ENNHK(3,fl-3)  -6D0*FNN  MK  (3  , ft- 2)  ♦ 12D0  *E  N N HK  (3,  M- 1 ) 

G -10DT*ENNHK(3,fl) ♦ 3 CO* ENNMK  (3, H*  1) ) /BON  3 

AJJH(2)  = (-AJJHK(2,K-3)  *6D0*AJJHK (2, K-2)  -1  8D0*  A J JMK  (2,  H-  1) 
fi  *1CD0*AJJHK(2,H)*3C0*AJJHK(2,H*1))/BrN1 

A JJ  H (1)  = (-AJJHK  ( 3,  r - 3 ) *4D0*AJJKK  ( 3,  H-2)  *6D0  *A  J J HK  ( 3 , H - 1 ) 

& -20D0*A  J JHK  ( 3,H)  ♦ 1 1DO*  A JJHK  ( 3,  F*  1) ) /BDf‘2 

GO  TO  85 

902  CONTINUE 

SNNH(2)  = (1lr0*ENNKK  (2, fl-4)  -56 DO*ENN HK ( 2, H - 3) ♦ 1 1 4D0*EN NHK  (2, H- 2) 
E -104D0*ENNFK  (2»  K- 1 ) *35D0*ENNNK  (2 , F)  )/Bl)N2 
FNNH(3)  = (3D0*ENNF.K  { 3,  fl-4)  -14D0*ENNFK  (3,  H-  3)  *24D0*ENNHK  (3, H-2) 

E -18D0*FNNHK  (3,  FI  — 1 ) *5D0*ENNNK  (3,  H)  ) / BDK  3 
AJJH  (2)  = ( 3 D 0 * A J J H K (2, R-4)  -1fcD0*AJJPK  (2, H-3)  *36D0*AJJHK(2, H-2) 

E -48D0*AJJHK  (2,H-1)  +25  DO  * A J J HK  (2, H) ) /3DN1 

AJJH  ( 3)  = ( 1 1D0*AJJ«K  \ 3, F-4)  -56P0*AJJ,,K  ( 3,  H-  3)  ♦ 11  4D0*AJJHK  (3,  H-2) 
C -104DO* AJJHK  (3, M- 1) ♦ 35D0* AJJHK  (3, H) ) /BDN2 
GO  TO  85 

801  CONTINUE 
T = H*9 

■ = ♦ 1 

IF (T.GT.TTU .ON.H.GT. HLIHK ) GO  TO  903 
DO  8 8«  J = 1 , 5 
888  DET  (J)  = ZE°0 
DO  811  J = 1 , 3 
ENE  (J) =ZEPO 

811  AJK  (J)  =ZEPO 
1 = 1 

812  CONTINUE 

Xfl  = (T-TTP) /TTQ 
TLU  = (T- TTL)  / (TTIJ-T) 

GII=GI  (I) 

SLU*.2D0/DSQRT (IDO* HI  (I) *TLO) 

yII=XSN*(XH* (1D0-XH) *GII*irO) /2D0 

yi2=Qll*cn 

y°I=1DG-5C0*QI2 

ySI=1DO-3DO*QI2 

OCI=yII * (2C0*0RI) 

IH=3 

I P (T . GT • TTUB2)  I H=4 
IP  (T.GT . TTU E ) IH=5 
IL=  1 

802  I Y= I YL- 3 

TL  = T* (IL-IH) *B 

803  IY=IY*1 

TB  = QK  (I Y)  *53(11)  *QI  I 
IF(TL.GT.TB)  GO  TO  803 
IP (IY.GE. IYO)  I Y=I Y V 
TB  = QK(I  Y)  *SS  (II)  *QI I 
IB  A = II- 2 
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XHB  = IY-  1 
X W C = I Y * 1 

TWA=0K (lb  A)  *SS  (IWA) 

T H b = O K (IWB)  ♦ SB  (IbB)  *0IX 

TUC  = UK (1WC)  ♦ SS  (TWC)  *yix 

Ta::=  (-THA  + ^no*  (ti«3*tw)  -■"« c)  /16do 

TWC= (TWA-27P0*  (TWB-TW)  -TWC) /24C0 

1HU=  (TW  A-TbB-TW*TWC)  /4P0 

T b!i  = (TL-TWH)/THF 

TUI= (-TWA»3DC*  (TWR-TW  ) *TWC)  /BOO 
TWJ=?MG/r«r 

rwK=2no*vwj*T«j-Twi/"wr 

AP  = . SD(j  ♦ Tbh*  (1DO-TWH*  ( TW,7  - T WK  *T  WH)  ) 

AP*'2  = AP-2D0 
AFr.!1=AP*AP-inO 
AL1=-AP*(AP-1DO)  *AP«2/6nr> 

AL2  = 0.5D0*AF  GP*Arf’2 

AL3=-0.  87)0*  A P*  (AP+1D0)  *A?!".2 

^LU  = AP*  Ar'sn/btO 

Dn  b04  ,7=1,  u 

JY= IY-J* J 

Z=  { ,1 Y - 1 ) *!I*Q 

SINZ=DSIN  (Z) 

C03Z=Dcr.‘j  (Z) 

B U ,7  Y = BU  (.IY) 

DE ( IL.J) = IDC-BDGY* ( X*COSZ-0II *51 NZ) 

AK1  (IL, J)  = (X*SINZ*OII  *COSZ)  *BIJY 

»K2  (IL.J)  = (FF  (8,  JY)  *0II*FF  (0,.1Y)  *0SI  *FF  (10,  JY)  ) *PtIJY 
A K 3 (I  L,  .1)  = (FF  (11,  JY)  ♦ OH*  FF  (1  2 , ,1 Y)  ♦0PI*FF  (1  3,  ,1Y ) ♦ 0C I*  FF  (14, JY)  ) 

6 *BUJY 

ES2  (XL,  J)  =FF<1,,IY)  ♦ QI I *FF  (2,JY)  *0-81  *FF  ( 3,  .7  Y) 

804  '‘Hi  ( I L , J ) = F F ( 4 , ,1 Y ) ♦ QII*F',’(S,JY)  ♦'.)?!*  FF  ( b,  .1 Y ) ♦ QC  T * FF  (7  , J Y ) 

DDA  (I’.)  =AL1*DL  (IL,  1 ) ♦ AI.  2 *DE  ( XL,  2)  ♦ AL3*D  F ( IL , ))  ♦AX4*Dc.(:L,4) 

DD AI =DDA (IL) 

ENA  (1 ,1 L)  =1 CO/DPAI 

E*7A  (2,TL)  =(AL1*E82(IL,  1l  t A L2* EH2  ( IL , 2)  fAL3*E«2(IL,3)  ♦ AL4  * ;;r  2 ( I L , P ) 
r.  l/CDAI 

E”A  ( 1,  I L)  = (M.1*’?"3(Il  , 1)  ♦ A L2* " “ 3 (II.  ,2)  ♦ A L 3*  EF  1 (IE,  1)  ♦ A 14  * ME  3 (II  , u ) 
r.  ) /P  E A I 

A. 7 A ( 1 , I L ) = ( A L 1 * A K 1 (II  , 1)  ♦ A L 2*  A K 1 (IL,  2)  ♦ AL  3 ♦ AK  1 (IL,  1)  *A14*AK1  ( I I.  , 4 ) 

f,  )/LllAr 

A,JA  (2,1  l)  = ( f LI*  AK2  ( IL,  1)  * AL2*  AK2  (IL,  w)  ♦ AL  3*  AK2  ( IT  ,3)  ♦ AL4*AK2  (IL,4) 
f,  ) /D  D A 1 

A.7A  (1,11  ) = ( A L 1 * A K 3 (IL,  1)  ♦ AL2*AK3  (IL,  2)  ♦ AL  3*AK3  (IL,  3)  ♦ AL4*AK3  (IL,4) 
C ) /DDA I 

IL= XL*  1 

IF (IL.LF. 8 ) GO  *0  802 
or  80S  ,1=1,  1 

FNK (J)  = ENS ( J) ♦SLU*FN A (J, IK) 

805  A.1K  ( ,1 ) = A ,1  K (.7)  ♦BLU*AJA  (0,18) 

IP  (Il.Xr  . ■*)  GO  TC  8 0b 

DEC (1)=  ( - 1 N A (2, 1) ♦ 1 SCO* (ENA (2,2)  *FN  A (2, 4) ) -J?CO*FNA (2,3) 

6 -El.  A (2,  S)) /3D«2 

DEO  (2)=  (-FNA  (3,  1)  *2D0*  (FN  A (3,  2)  - ENA  ( 3,4))  +5NA  (3,5) ) /BDN3 
DE0  (3) = ( A.7  A (2,  1)  -8  DO*  (A/I  A (2 , 2)  - A,1A(2, 4 ) ) - A,1  A (2  , S)  ) /B  DN  1 
CEQ («) = (-AJ A (3,  1) ♦ 1 6D0* (A J A (3,2) +AJA  (3,4) ) - 30D0*AJA (1, 3) 

C - AJ  A ( 3 , S) ) /BON 2 

r,n  jr  flOB 

00b  IF(IM.E0.5)  go  to  807 

DEO  (1)  = (-LNA  (2,  1)  *4  DO*  EM  A (2,2)  tbDO*E*’A  (2,  3)  -20D0*ENA  (2,4) 

F,  *11  DO*FNA (2,5)  ) /BDS2 
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DEO (2) = (ENA (3, 1) -6D0*ENA (3, 2) +12D0*ENA (3, 3) -1 ODO*ENA (3, 4) 
e ♦3l)0*ENA  (3,5)  ) / BDN3 

DEO  (3)  = (- AJA<2, 1)  *-6P0*AJA  (2,2) -18D0*AJA  (2, 3) ♦10D0*AJA  (2,4) 

6 ♦ 3DO*AJ A (2, 5) ) /BDN  1 

DE<J  (4)  = (-AJA  (3,  1)  ♦ 4 DO  * AJA  (3,2)  ♦6DO*AJA  (3,  3)  -20D0*AJA  (3,4) 

6 ♦ 1 1 DO*AJA  (3, 5) ) /BDN2 

GO  TO  S308 
BC7  CONTINUE 

DEQ  (1)  = (11D">*ENA  (2,  1)  -5bP  0*  E”  A (2 , 2)  ♦ 1 14  DO* ENA  (2 , 3)  - 1 04D0*  BN  A (2, 4) 
6 ♦ )5D0*FNA(2,5)  ) / BCN2 

DE2  (2)  = ( 3D'1*FNA  (3,  1)  -14D0*ENA  (3,  2)  ♦ 24  DO  *FNA  ( 3 , 3)  - 1 8D0*EN  A (3,4) 
6 + 5D0*E>'A  (3,8))  /BDN3 

DEO  (3) = ( 3D0*AJA  (2, 1)  -16D0* AJ A (2, 2)  ♦ 36 DO* A JA ( 2 , 3) - 48D0*AJA (2 , 4) 
f,  ♦25D0*AJA  (2, 5)  ) /BDN1 

DEO  (4)  = (3 IDO *AJ A (3,  1)  -56B0* AJ A(3,2)*114D0*AJA().3)-,C4D0*AJA(3,4) 
& ♦ 35  DC  * A J A (3,5))  /BDN2 

808  C 0 I N U E 

DO  804  J=  1 , 4 

809  DEB  (J)  =DEP  (J)  ♦ SLtl*DEQ  (J) 

r = i ♦ i 

I F ( I . LE  . 5)  GO  T o «12 

DO  810  J=1, 3 

ENNMK  (J  , ■*)  = '‘ACQ-  *EN  K ( J) 

810  A J J UK (J,N) = 0. 90Ub7D0*PACOM*AJK  (0) 

ENNfl  (1)  =ENN  ■'K  (1  ,*) 

ENNf!  (2)  =DEB  (1)  ♦FACC" 

ENNM  (3)  = DEI  (2)  *FACQH 

ajjm  (i)  =AJjrK  ( 1 , r.) 

AJJH(2)  =0.9C467D0*DES  (3)  *FACOM 
A J J R ( 3)  =0.9C.4b7DC*DER  (4)  *FACy* 

85  CONTINUE 

SU-N  = ENN*  (1)  ♦ FNV*  (?)  *ENN"  (3) 

FNN  (fl) =?N  N ( fi)  ♦ .5 DO*  (SUI*N*ENF  ( E)  ) 

ENP  (.1)  = SURN 

IF  (YY. EC. 3010)  FNP(«)  = SUM N ♦ SU 1 N 

Stl-J  = A JJN  (1)  ♦AJJW  (2)  ♦AJJP.  (3) 

AJJ  (M)=AJJ(Pi)  f .5  DO*  (SUMJ+AJF (1)  ) 

A.1  P (N)  = SUNJ 

IF  l YY.EQ.  30 DO)  AJP(fl)  = SURJ*SUflJ 
IF((2*((R-1)/2)  . FO. F-1 . AN  P. T. LE.2.01D0)  .0s. 

C (5*((!4-1)/5)  .EQ ."-1 .AKD.T.GT.2.01D0) ) 

6PRINT  5 52.T,  (AJJR(I) ,1=1,3)  ,SUMJ, AJJ (tl) 

6 , (ENKH(I) ,1=1, 3) ,SDHN,ENN(H) 

552  FOPRAT (F8.2,2  (1 X,4F9. 5,  F10.5)) 

IF (N. EQ . 1 ) GO  TO  91 
IF (R.GF . MLX)  GO  TO  801 
IF (l.tT.FLIFK)  GO  TO  92 

903  CONTINUE 

IP  (PUX.  Ey.RLIMK)  GO  TO  905 

904  CONTINUE 

TVU  = .5D0*  (TTV-T) / (TTV-TTU) 

ENN  (M)=ENN(F)  ♦ TVU*ENP(R)  * FPC  3 
AJJ  (N) = AJJ ( 8)  ♦TVU*AJP  (H) 

IF  ( (2*(  (F-1)/2) .EQ.R-1.AND.T. LE.2.01D0)  .OR. 

6 (5*  ( (1*-1) /5)  . Ey  .11-1. AND. T .GT.  2. 01  DO)) 

&PFINT  855, T, AJJ ( fl)  , ENN  (N) 

855  PORHAT(F8.2,2(37X,  F10.5)) 
t=r*b 
b = r*i 

IF(H.LE.nVX)  GO  TO  904 


n n 
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905  T?V 

nvx 

F. 

ENT 


S UBKn!IT  I S'  L PAFT7 
*******  PA'-'T  7 rEVIS’:[; 

******  U-POI  "T  V OF SIGN 

IMPLICIT  F S A L * 8 (A-H,K,0-Z) 

common  (c,o,  121 ) , TiSN  ( i ii)  , ajj  (i  1 1) , rr  (5,  ui) 

6,  F (200)  ,CC<  121)  ,SS(121)  ,T1  (5,  121)  ,DP(5,  121)  ,*0(50) 
fi.FNNMK  (3,1  1 1)  , A.UKK  ( 3 , 1 1 1 ) , . P i ( 5,  1 2 1)  , FN2  ( 8,  1 2 1 ) , EN  3 ( 5,  1 2 1 ) 

6,  ENN*  (*)  , A.MK  (3)  , AJ  1 (5,12  1)  , A.J2  (5, 1 21)  , A.J3  (5,  12  1) 

6,  F-F(3)  ,AJK(M  (1,9),  AJA  (1,9) 

6 , D 0 A (9)  , 01  («)  ,0C  (J)  ,0°(9)  ,0S(°)  ,1<U  (121)  ,0E0(5)  , Ot'F  (c) 

6, OK  (121)  , FT  (14. 121)  ,51  (5)  , HI  ( 5)  ,DE(5,4)  ,AK1  (5,4) 

6, AK2  (5,4)  ,AK3(5,4) , FM2  (5, 4)  ,r *J  (5,4) , ENM(3)  ,"X2  (4) , EK 3 (4) 

6,  Evo  (ii  i)  , djp  (i  i i) 

C OF  MOV  HH,DLLH,YY,YYPEL,VV,r,W,G6,KK,PB,E’'Y,X,XSN,P.Hf:,B,C,YO 
6 , T rP,TTQ,TTL,i’T'J,T’'V 

6,  "’LIK.MLIPK.IH.LV,!  YML  ,TYU,TYL,  LGO,  ] YV,  «VX 
3(ZDUM)  =33*rFXI  (-KK*ZDUM) 

IMI  ( ZD'J*’ ) = ( (('.  (ZP'JF)  - 1D0)  * ( KG  ♦ 1 DO)  / ( (G  (ZOO*;)  ♦ IDO)  * (GG-1  DO)  ) ) 

1 **2  . 1 64  DC 

A (ZOOM)  =2 . 1<i«D0*KK*G  (ZFUf.)  / (3  (ZP'tM)  **2-  IDO) 

0 (ZD'J*)  =5. 885690  C*D  SOFT  (G  (ZDDV)  **2-1  DO)  *!M1  (ZDUK)  /3B 

DATA  ZrrO,r,C/0D0, 3. 1 4155265353979D0, 1 . 74532925149U j JO- 2/ 

6 , v*  , Qv/0  • 95 1 056  5 16  2 95  1 5D0  , C • 5 £778  525229  24  7D  0/ 

700  CONTINUE 

DO  701  1=  1,20  0 
7C1  f(I)  =ZEF<~ 

DO  702  M Y = 1 , 1 2 1 
Dn  702  1=1,50 
702  "I  (I, MY)  =ZFrC 
IV2=I.V*!.V 
LV  3 = LV2  ♦ L V 
j = v n 
XGO  = 1 
X Y = 1 

IF  (IYKL.FO. 1) 

6 SETUP'! 

3G  TO  ( 6 0 2 , 3 2 , o 4 0 ) ,L3C 
3 0 c n '■  T I 'IU  0 

1 = 1 

32  C(I)  = F(I*LV3) 

1 = 1*1 

IF  (I.LE.T-V)  GO  TO  3 2 

33  Y = il  ♦ Y 
IO''  = 2 

X Y = I Y ♦ 1 

GO  (602,  *>32,640)  ,LGC 
35  CONTINUE 
1=1 

37  F (I  *LV)  =F  ( I ♦ LV 3) 

1 = 1*1 

IF  (I.LE.LV)  30  TO  3 7 
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IF  (IYNL.GT.  2)  GO  TO  701 
1 = 1 

41  El  (1 , 2)  =H*Q*  (?  (I  J *F  (I*LV)  )*.5D0 
1 = 1*1 

IF  (I.LE.LV)  GO  TO  4 1 
FETU3N 

703  Y = Y*H 
IGO=l 

I Y= I Y* 1 

GO  TO  (602,632,640)  ,LGC 
40  CONTINUE 
1=1 

42  F (I ♦ LV2 ) = F ( I ♦ t>  V 3 > 

1 = 1*1 

IF  (I.LE.LV)  GO  TO  42 
I P ( IYML . GT. 3)  GO  T"  704 
1 = 1 

4 4 El  (1,2)  =H*0*  (5D0*F(I) *8D0*F  (I*LV) -F (I*LV2) ) /I  2 DO 
El  (1,3)  =U*0*(F(I)*4D0*F(I*LV)  *F  (I  + LV2)  ) /3D0 
1 = 1*1 

IF  (I.LE.LV)  GO  TO  44 
r.ETUPN 

704  Y=Y*H 
I GO  = 4 

I Y = I Y* 1 

GO  TO  (602,632,640) ,LGC 
45  COMTINUE 
1 = 1 

47  El  (1,2  ) =H*0*  (9D0*F  (I) *19D0*F  (I*LV) -5DC*F  (I*LV2)*F(I*Ljr3))/24D0 
1=1*1 

IF  (I.LE.LV)  GO  TO  47 
50  CONTINUE 
1=1 

4 8 El (I,IY-1) =EI  (I, IV- 2)  *H*0*  (-F  (I)  ♦ 13D0* (F(I*LV) *F (I*LV2) ) 

6 -F(I*LV3))/24D0 

1 = 1*1 

IF  (I.LE.LV)  GO  TO  48 
1=1 

IF ( Y*H. GT. W)  GO  TO  60 
52  F (I) =F( I*L V) 

F (I *LV)  =F  (I *LV2) 

F (I +LV2) =F(I*LV3) 

1 = 1*1 

IF(I.LE.LV)  GO  TO  52 
Y = Y*H 
IGO  = 5 
I Y= I Y*  1 

GO  TO  (602,632,640) ,LGO 
55  GO  TO  50 
60  CONTINUE 

62  El  (1,1  =EI(I,IY-1)  *H*C*(F(I)  -5D0*F(I*LV)  ♦ 1 9D0* F (I* LV 2) 

6 *9>,  - P (T*LV3)  )/24D0 

1 = 1*1 

IF(I.LE.LV)  GO  TO  62 
FETUEN 

C******* 

602  CONTINUE 
1 = 1 
Z = Y*0 
OOZ*0(Z> 


Ft  FT  8 
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f (i ♦ *>)  = ooz*rcos  (Z) 

F ( 1 ♦ *7 ) = OCZ*DSIW(Z) 

00  TO  (30, 35, 40, 45, 55)  , I.JO 

(-•••••**  rftPT  38 

6 32  CONTINUE 
Z = Y *0 
AOZ  = A (Z) 

II UO Z = UH  (Z) 

IIUZS  = 'IUOZ*UU('Z 
00(1)  = AOZ/(tlUZS*IJUOZ) 

00(22)  = Aoz*UUZtl 

1=1 

continue 

ooi=QO(i) 

CCI  = CC  (IY) 

SSI=S3  (I Y) 

F c I ♦ 1 26 ) =U0T*PP(1,IY) 

F (1*127) =F  (1*126) *CCI 
F ( I ♦ 1 28  ) = F ( x ♦ 1 2 6 ) *SSI 
F (1*12°) =1  ( I*  1 27) *3  SI 
F (1*1  30) =F  ( 1*127) *CC1 
F(I*131)  =F(I*12fl)  *5Si. 

F (1*132) =Cv  I*Pr  ( 2 , 1 Y ) 

F (1*1 33) =F  ( T* 1 32) *CCi 
F (I*  1 34) =F  ( I* 1 32) *SST 
r (I*135)  = F (I*13J)*SSI 
F (I *136) =F  (1*133) *CC1 
F (I *1  37) =F  (1*134) *SSI 
F (1*138) =wj:*pr  ( 3, 1 Y) 

F (I *139) =F( I* 1 38) *CCI 
F (I ♦ 14C) =F  ( I* 1 38) *SSI 
f (1*141)  = 00  1*0P  (4,1 Y) 
f (1*142) =F  (1*141) *CCI 
F(I*143)=F(:*141)*5S1 
F (1*144)  =0UI»PP ( 6 , I Y ) 

F (I *145) = F (1*144) *CCI 
F(X*146)  =F  ( 1*144) * 5 S I 
1=1*21 

IF  (I.  IT.  2 3)  00  rr  3 8 
00  TO  (30,35,40,45,55) ,100 

C*******  PAPT  46 

640  CONTINUE 
Z = Y*0 
APZ  = A (Z) 

UU0Z=UI1  (Z) 

UU7.S  = UU0Z*1JU0Z 

uuzc=uuoz*inizs 

00  (1) = AOZ*UIIOZ 
00(11) = AOZ/UUZS 
00(21)  = AOZ*UIIZC 

00  (31)  = AOZ/UUZC 
OQ  (41) =AOZ*UUZS 

1 = 1 

46  CONTINUE 
001=00(1) 

SSI =SS ( I Y) 

CCI =CC ( IY ) 

F (I ♦ 150) *0QI 
F (1*151 ) =QQI 


•CCI 


F (I ♦ 1 52 ) =QQ  X ♦SSI 
F (I* 153) =P(I*151) *SSI 
F (1*154) = F (1*151) *CCI 
F ( I ♦ 155) =F(I*152) *SSI 
F (1*156) =F (1*154) *CCI 
F (1*157) =F  (1*155) ♦CCI 
F (1*158) =F  ( I* 1 54) *SSI 
F (1*159) = F ( I ♦ 1 55) *SSI 
1=1*10 

IF(I.LT.42)  GO  TO  46 

GO  TO  (30,35,40,45,55)  ,IGC 

END 
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